/*
 * Copyright 2014 Actelion Pharmaceuticals Ltd., Gewerbestrasse 16, CH-4123 Allschwil, Switzerland
 *
 * This file is part of DataWarrior.
 * 
 * DataWarrior is free software: you can redistribute it and/or modify it under the terms of the
 * GNU General Public License as published by the Free Software Foundation, either version 3 of
 * the License, or (at your option) any later version.
 * 
 * DataWarrior is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
 * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 * You should have received a copy of the GNU General Public License along with DataWarrior.
 * If not, see http://www.gnu.org/licenses/.
 *
 * @author Joel Freyss
 */
package com.actelion.research.forcefield;


/**
 * Derived from apache math package
 *
 */
public class FastMath {
	private static final double SAFE_MIN = 0x1.0p-1022;

	/** Archimede's constant PI, ratio of circle circumference to diameter. */
	public static final double PI = 105414357.0 / 33554432.0 + 1.984187159361080883e-9;

	/** Napier's constant e, base of the natural logarithm. */
	public static final double E = 2850325.0 / 1048576.0 + 8.254840070411028747e-8;

	/**
	 * Exponential evaluated at integer values, exp(x) = expIntTableA[x + 750] +
	 * expIntTableB[x+750].
	 */
	private static final double EXP_INT_TABLE_A[] = new double[1500];

	/**
	 * Exponential evaluated at integer values, exp(x) = expIntTableA[x + 750] +
	 * expIntTableB[x+750]
	 */
	private static final double EXP_INT_TABLE_B[] = new double[1500];

	/**
	 * Exponential over the range of 0 - 1 in increments of 2^-10 exp(x/1024) =
	 * expFracTableA[x] + expFracTableB[x].
	 */
	private static final double EXP_FRAC_TABLE_A[] = new double[1025];

	/**
	 * Exponential over the range of 0 - 1 in increments of 2^-10 exp(x/1024) =
	 * expFracTableA[x] + expFracTableB[x].
	 */
	private static final double EXP_FRAC_TABLE_B[] = new double[1025];

	/** Factorial table, for Taylor series expansions. */
	private static final double FACT[] = new double[20];

	/**
	 * Extended precision logarithm table over the range 1 - 2 in increments of
	 * 2^-10.
	 */
	private static final double LN_MANT[][] = new double[1024][];

	/** log(2) (high bits). */
	private static final double LN_2_A = 0.693147063255310059;

	/** log(2) (low bits). */
	private static final double LN_2_B = 1.17304635250823482e-7;

	/** Coefficients for slowLog. */
	private static final double LN_SPLIT_COEF[][] = { { 2.0, 0.0 }, { 0.6666666269302368, 3.9736429850260626E-8 }, { 0.3999999761581421, 2.3841857910019882E-8 }, { 0.2857142686843872, 1.7029898543501842E-8 },
			{ 0.2222222089767456, 1.3245471311735498E-8 }, { 0.1818181574344635, 2.4384203044354907E-8 }, { 0.1538461446762085, 9.140260083262505E-9 }, { 0.13333332538604736, 9.220590270857665E-9 }, { 0.11764700710773468, 1.2393345855018391E-8 },
			{ 0.10526403784751892, 8.251545029714408E-9 }, { 0.0952233225107193, 1.2675934823758863E-8 }, { 0.08713622391223907, 1.1430250008909141E-8 }, { 0.07842259109020233, 2.404307984052299E-9 }, { 0.08371849358081818, 1.176342548272881E-8 },
			{ 0.030589580535888672, 1.2958646899018938E-9 }, { 0.14982303977012634, 1.225743062930824E-8 }, };

	/** Coefficients for log, when input 0.99 < x < 1.01. */
	private static final double LN_QUICK_COEF[][] = { { 1.0, 5.669184079525E-24 }, { -0.25, -0.25 }, { 0.3333333134651184, 1.986821492305628E-8 }, { -0.25, -6.663542893624021E-14 }, { 0.19999998807907104, 1.1921056801463227E-8 },
			{ -0.1666666567325592, -7.800414592973399E-9 }, { 0.1428571343421936, 5.650007086920087E-9 }, { -0.12502530217170715, -7.44321345601866E-11 }, { 0.11113807559013367, 9.219544613762692E-9 }, };

	/** Coefficients for log in the range of 1.0 < x < 1.0 + 2^-10. */
	private static final double LN_HI_PREC_COEF[][] = { { 1.0, -6.032174644509064E-23 }, { -0.25, -0.25 }, { 0.3333333134651184, 1.9868161777724352E-8 }, { -0.2499999701976776, -2.957007209750105E-8 },
			{ 0.19999954104423523, 1.5830993332061267E-10 }, { -0.16624879837036133, -2.6033824355191673E-8 } };

	/** Sine table (high bits). */
	private static final double SINE_TABLE_A[] = new double[14];

	/** Sine table (low bits). */
	private static final double SINE_TABLE_B[] = new double[14];

	/** Cosine table (high bits). */
	private static final double COSINE_TABLE_A[] = new double[14];

	/** Cosine table (low bits). */
	private static final double COSINE_TABLE_B[] = new double[14];

	/** Tangent table, used by atan() (high bits). */
	private static final double TANGENT_TABLE_A[] = new double[14];

	/** Tangent table, used by atan() (low bits). */
	private static final double TANGENT_TABLE_B[] = new double[14];

	/** Bits of 1/(2*pi), need for reducePayneHanek(). */
	private static final long RECIP_2PI[] = new long[] { (0x28be60dbL << 32) | 0x9391054aL, (0x7f09d5f4L << 32) | 0x7d4d3770L, (0x36d8a566L << 32) | 0x4f10e410L, (0x7f9458eaL << 32) | 0xf7aef158L, (0x6dc91b8eL << 32) | 0x909374b8L,
			(0x01924bbaL << 32) | 0x82746487L, (0x3f877ac7L << 32) | 0x2c4a69cfL, (0xba208d7dL << 32) | 0x4baed121L, (0x3a671c09L << 32) | 0xad17df90L, (0x4e64758eL << 32) | 0x60d4ce7dL, (0x272117e2L << 32) | 0xef7e4a0eL,
			(0xc7fe25ffL << 32) | 0xf7816603L, (0xfbcbc462L << 32) | 0xd6829b47L, (0xdb4d9fb3L << 32) | 0xc9f2c26dL, (0xd3d18fd9L << 32) | 0xa797fa8bL, (0x5d49eeb1L << 32) | 0xfaf97c5eL, (0xcf41ce7dL << 32) | 0xe294a4baL, 0x9afed7ecL << 32 };

	/** Bits of pi/4, need for reducePayneHanek(). */
	private static final long PI_O_4_BITS[] = new long[] { (0xc90fdaa2L << 32) | 0x2168c234L, (0xc4c6628bL << 32) | 0x80dc1cd1L };

	/**
	 * Eighths. This is used by sinQ, because its faster to do a table lookup
	 * than a multiply in this time-critical routine
	 */
	private static final double EIGHTHS[] = { 0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625 };

	/** Table of 2^((n+2)/3) */
	private static final double CBRTTWO[] = { 0.6299605249474366, 0.7937005259840998, 1.0, 1.2599210498948732, 1.5874010519681994 };

	/*
	 * There are 52 bits in the mantissa of a double. For additional precision,
	 * the code splits double numbers into two parts, by clearing the low order
	 * 30 bits if possible, and then performs the arithmetic on each half
	 * separately.
	 */

	/**
	 * 0x40000000 - used to split a double into two parts, both with the low
	 * order bits cleared. Equivalent to 2^30.
	 */
	private static final long HEX_40000000 = 0x40000000L; // 1073741824L

	/** Mask used to clear low order 30 bits */
	private static final long MASK_30BITS = -1L - (HEX_40000000 - 1); // 0xFFFFFFFFC0000000L;

	/**
	 * 2^52 - double numbers this large must be integral (no fraction) or NaN or
	 * Infinite
	 */
	private static final double TWO_POWER_52 = 4503599627370496.0;

	// Initialize tables
	static {
		int i;

		// Generate an array of factorials
		FACT[0] = 1.0;
		for (i = 1; i < FACT.length; i++) {
			FACT[i] = FACT[i - 1] * i;
		}

		double tmp[] = new double[2];
		double recip[] = new double[2];

		// Populate expIntTable
		for (i = 0; i < 750; i++) {
			expint(i, tmp);
			EXP_INT_TABLE_A[i + 750] = tmp[0];
			EXP_INT_TABLE_B[i + 750] = tmp[1];

			if (i != 0) {
				// Negative integer powers
				splitReciprocal(tmp, recip);
				EXP_INT_TABLE_A[750 - i] = recip[0];
				EXP_INT_TABLE_B[750 - i] = recip[1];
			}
		}

		// Populate expFracTable
		for (i = 0; i < EXP_FRAC_TABLE_A.length; i++) {
			slowexp(i / 1024.0, tmp);
			EXP_FRAC_TABLE_A[i] = tmp[0];
			EXP_FRAC_TABLE_B[i] = tmp[1];
		}

		// Populate lnMant table
		for (i = 0; i < LN_MANT.length; i++) {
			double d = Double.longBitsToDouble((((long) i) << 42) | 0x3ff0000000000000L);
			LN_MANT[i] = slowLog(d);
		}

		// Build the sine and cosine tables
		buildSinCosTables();
	}

	/**
	 * Private Constructor
	 */
	private FastMath() {
	}

	// Generic helper methods

	/**
	 * Get the high order bits from the mantissa. Equivalent to adding and
	 * subtracting HEX_40000 but also works for very large numbers
	 *
	 * @param d
	 *            the value to split
	 * @return the high order part of the mantissa
	 */
	private static double doubleHighPart(double d) {
		if (d > -SAFE_MIN && d < SAFE_MIN) {
			return d; // These are un-normalised - don't try to convert
		}
		long xl = Double.doubleToLongBits(d);
		xl = xl & MASK_30BITS; // Drop low order bits
		return Double.longBitsToDouble(xl);
	}

	/**
	 * Compute the square root of a number.
	 * <p>
	 * <b>Note:</b> this implementation currently delegates to {@link Math#sqrt}
	 * 
	 * @param a
	 *            number on which evaluation is done
	 * @return square root of a
	 */
	public static double sqrt(final double a) {
		return a<=0? a: Math.sqrt(a);
	}

	/**
	 * Compute the hyperbolic cosine of a number.
	 * 
	 * @param x
	 *            number on which evaluation is done
	 * @return hyperbolic cosine of x
	 */
	public static double cosh(double x) {
		if (x != x) {
			return x;
		}

		if (x > 20.0) {
			return exp(x) / 2.0;
		}

		if (x < -20) {
			return exp(-x) / 2.0;
		}

		double hiPrec[] = new double[2];
		if (x < 0.0) {
			x = -x;
		}
		exp(x, 0.0, hiPrec);

		double ya = hiPrec[0] + hiPrec[1];
		double yb = -(ya - hiPrec[0] - hiPrec[1]);

		double temp = ya * HEX_40000000;
		double yaa = ya + temp - temp;
		double yab = ya - yaa;

		// recip = 1/y
		double recip = 1.0 / ya;
		temp = recip * HEX_40000000;
		double recipa = recip + temp - temp;
		double recipb = recip - recipa;

		// Correct for rounding in division
		recipb += (1.0 - yaa * recipa - yaa * recipb - yab * recipa - yab * recipb) * recip;
		// Account for yb
		recipb += -yb * recip * recip;

		// y = y + 1/y
		temp = ya + recipa;
		yb += -(temp - ya - recipa);
		ya = temp;
		temp = ya + recipb;
		yb += -(temp - ya - recipb);
		ya = temp;

		double result = ya + yb;
		result *= 0.5;
		return result;
	}

	/**
	 * Compute the hyperbolic sine of a number.
	 * 
	 * @param x
	 *            number on which evaluation is done
	 * @return hyperbolic sine of x
	 */
	public static double sinh(double x) {
		boolean negate = false;
		if (x != x) {
			return x;
		}

		if (x > 20.0) {
			return exp(x) / 2.0;
		}

		if (x < -20) {
			return -exp(-x) / 2.0;
		}

		if (x == 0) {
			return x;
		}

		if (x < 0.0) {
			x = -x;
			negate = true;
		}

		double result;

		if (x > 0.25) {
			double hiPrec[] = new double[2];
			exp(x, 0.0, hiPrec);

			double ya = hiPrec[0] + hiPrec[1];
			double yb = -(ya - hiPrec[0] - hiPrec[1]);

			double temp = ya * HEX_40000000;
			double yaa = ya + temp - temp;
			double yab = ya - yaa;

			// recip = 1/y
			double recip = 1.0 / ya;
			temp = recip * HEX_40000000;
			double recipa = recip + temp - temp;
			double recipb = recip - recipa;

			// Correct for rounding in division
			recipb += (1.0 - yaa * recipa - yaa * recipb - yab * recipa - yab * recipb) * recip;
			// Account for yb
			recipb += -yb * recip * recip;

			recipa = -recipa;
			recipb = -recipb;

			// y = y + 1/y
			temp = ya + recipa;
			yb += -(temp - ya - recipa);
			ya = temp;
			temp = ya + recipb;
			yb += -(temp - ya - recipb);
			ya = temp;

			result = ya + yb;
			result *= 0.5;
		} else {
			double hiPrec[] = new double[2];
			expm1(x, hiPrec);

			double ya = hiPrec[0] + hiPrec[1];
			double yb = -(ya - hiPrec[0] - hiPrec[1]);

			/* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
			double denom = 1.0 + ya;
			double denomr = 1.0 / denom;
			double denomb = -(denom - 1.0 - ya) + yb;
			double ratio = ya * denomr;
			double temp = ratio * HEX_40000000;
			double ra = ratio + temp - temp;
			double rb = ratio - ra;

			temp = denom * HEX_40000000;
			double za = denom + temp - temp;
			double zb = denom - za;

			rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;

			// Adjust for yb
			rb += yb * denomr; // numerator
			rb += -ya * denomb * denomr * denomr; // denominator

			// y = y - 1/y
			temp = ya + ra;
			yb += -(temp - ya - ra);
			ya = temp;
			temp = ya + rb;
			yb += -(temp - ya - rb);
			ya = temp;

			result = ya + yb;
			result *= 0.5;
		}

		if (negate) {
			result = -result;
		}

		return result;
	}

	/**
	 * Compute the hyperbolic tangent of a number.
	 * 
	 * @param x
	 *            number on which evaluation is done
	 * @return hyperbolic tangent of x
	 */
	public static double tanh(double x) {
		boolean negate = false;

		if (x != x) {
			return x;
		}

		if (x > 20.0) {
			return 1.0;
		}

		if (x < -20) {
			return -1.0;
		}

		if (x == 0) {
			return x;
		}

		if (x < 0.0) {
			x = -x;
			negate = true;
		}

		double result;
		if (x >= 0.5) {
			double hiPrec[] = new double[2];
			// tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
			exp(x * 2.0, 0.0, hiPrec);

			double ya = hiPrec[0] + hiPrec[1];
			double yb = -(ya - hiPrec[0] - hiPrec[1]);

			/* Numerator */
			double na = -1.0 + ya;
			double nb = -(na + 1.0 - ya);
			double temp = na + yb;
			nb += -(temp - na - yb);
			na = temp;

			/* Denominator */
			double da = 1.0 + ya;
			double db = -(da - 1.0 - ya);
			temp = da + yb;
			db += -(temp - da - yb);
			da = temp;

			temp = da * HEX_40000000;
			double daa = da + temp - temp;
			double dab = da - daa;

			// ratio = na/da
			double ratio = na / da;
			temp = ratio * HEX_40000000;
			double ratioa = ratio + temp - temp;
			double ratiob = ratio - ratioa;

			// Correct for rounding in division
			ratiob += (na - daa * ratioa - daa * ratiob - dab * ratioa - dab * ratiob) / da;

			// Account for nb
			ratiob += nb / da;
			// Account for db
			ratiob += -db * na / da / da;

			result = ratioa + ratiob;
		} else {
			double hiPrec[] = new double[2];
			// tanh(x) = expm1(2x) / (expm1(2x) + 2)
			expm1(x * 2.0, hiPrec);

			double ya = hiPrec[0] + hiPrec[1];
			double yb = -(ya - hiPrec[0] - hiPrec[1]);

			/* Numerator */
			double na = ya;
			double nb = yb;

			/* Denominator */
			double da = 2.0 + ya;
			double db = -(da - 2.0 - ya);
			double temp = da + yb;
			db += -(temp - da - yb);
			da = temp;

			temp = da * HEX_40000000;
			double daa = da + temp - temp;
			double dab = da - daa;

			// ratio = na/da
			double ratio = na / da;
			temp = ratio * HEX_40000000;
			double ratioa = ratio + temp - temp;
			double ratiob = ratio - ratioa;

			// Correct for rounding in division
			ratiob += (na - daa * ratioa - daa * ratiob - dab * ratioa - dab * ratiob) / da;

			// Account for nb
			ratiob += nb / da;
			// Account for db
			ratiob += -db * na / da / da;

			result = ratioa + ratiob;
		}

		if (negate) {
			result = -result;
		}

		return result;
	}

	/**
	 * Compute the inverse hyperbolic cosine of a number.
	 * 
	 * @param a
	 *            number on which evaluation is done
	 * @return inverse hyperbolic cosine of a
	 */
	public static double acosh(final double a) {
		return FastMath.log(a + FastMath.sqrt(a * a - 1));
	}

	/**
	 * Compute the inverse hyperbolic sine of a number.
	 * 
	 * @param a
	 *            number on which evaluation is done
	 * @return inverse hyperbolic sine of a
	 */
	public static double asinh(double a) {

		boolean negative = false;
		if (a < 0) {
			negative = true;
			a = -a;
		}

		double absAsinh;
		if (a > 0.167) {
			absAsinh = FastMath.log(FastMath.sqrt(a * a + 1) + a);
		} else {
			final double a2 = a * a;
			if (a > 0.097) {
				absAsinh = a
						* (1 - a2
								* (1 / 3.0 - a2
										* (1 / 5.0 - a2 * (1 / 7.0 - a2 * (1 / 9.0 - a2 * (1.0 / 11.0 - a2 * (1.0 / 13.0 - a2 * (1.0 / 15.0 - a2 * (1.0 / 17.0) * 15.0 / 16.0) * 13.0 / 14.0) * 11.0 / 12.0) * 9.0 / 10.0) * 7.0 / 8.0) * 5.0 / 6.0)
										* 3.0 / 4.0) / 2.0);
			} else if (a > 0.036) {
				absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0 - a2 * (1 / 7.0 - a2 * (1 / 9.0 - a2 * (1.0 / 11.0 - a2 * (1.0 / 13.0) * 11.0 / 12.0) * 9.0 / 10.0) * 7.0 / 8.0) * 5.0 / 6.0) * 3.0 / 4.0) / 2.0);
			} else if (a > 0.0036) {
				absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0 - a2 * (1 / 7.0 - a2 * (1 / 9.0) * 7.0 / 8.0) * 5.0 / 6.0) * 3.0 / 4.0) / 2.0);
			} else {
				absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0) * 3.0 / 4.0) / 2.0);
			}
		}

		return negative ? -absAsinh : absAsinh;

	}

	/**
	 * Compute the inverse hyperbolic tangent of a number.
	 * 
	 * @param a
	 *            number on which evaluation is done
	 * @return inverse hyperbolic tangent of a
	 */
	public static double atanh(double a) {

		boolean negative = false;
		if (a < 0) {
			negative = true;
			a = -a;
		}

		double absAtanh;
		if (a > 0.15) {
			absAtanh = 0.5 * FastMath.log((1 + a) / (1 - a));
		} else {
			final double a2 = a * a;
			if (a > 0.087) {
				absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0 + a2 * (1.0 / 7.0 + a2 * (1.0 / 9.0 + a2 * (1.0 / 11.0 + a2 * (1.0 / 13.0 + a2 * (1.0 / 15.0 + a2 * (1.0 / 17.0)))))))));
			} else if (a > 0.031) {
				absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0 + a2 * (1.0 / 7.0 + a2 * (1.0 / 9.0 + a2 * (1.0 / 11.0 + a2 * (1.0 / 13.0)))))));
			} else if (a > 0.003) {
				absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0 + a2 * (1.0 / 7.0 + a2 * (1.0 / 9.0)))));
			} else {
				absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0)));
			}
		}

		return negative ? -absAtanh : absAtanh;

	}

	/**
	 * Compute the signum of a number. The signum is -1 for negative numbers, +1
	 * for positive numbers and 0 otherwise
	 * 
	 * @param a
	 *            number on which evaluation is done
	 * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
	 */
	public static double signum(final double a) {
		return (a < 0.0) ? -1.0 : ((a > 0.0) ? 1.0 : a); // return +0.0/-0.0/NaN
															// depending on a
	}

	/**
	 * Compute the signum of a number. The signum is -1 for negative numbers, +1
	 * for positive numbers and 0 otherwise
	 * 
	 * @param a
	 *            number on which evaluation is done
	 * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
	 */
	public static float signum(final float a) {
		return (a < 0.0f) ? -1.0f : ((a > 0.0f) ? 1.0f : a); // return
																// +0.0/-0.0/NaN
																// depending on
																// a
	}

	/**
	 * Compute next number towards positive infinity.
	 * 
	 * @param a
	 *            number to which neighbor should be computed
	 * @return neighbor of a towards positive infinity
	 */
	public static double nextUp(final double a) {
		return nextAfter(a, Double.POSITIVE_INFINITY);
	}

	/**
	 * Compute next number towards positive infinity.
	 * 
	 * @param a
	 *            number to which neighbor should be computed
	 * @return neighbor of a towards positive infinity
	 */
	public static float nextUp(final float a) {
		return nextAfter(a, Float.POSITIVE_INFINITY);
	}

	/**
	 * Returns a pseudo-random number between 0.0 and 1.0.
	 * <p>
	 * <b>Note:</b> this implementation currently delegates to
	 * {@link Math#random}
	 * 
	 * @return a random number between 0.0 and 1.0
	 */
	public static double random() {
		return Math.random();
	}

	/**
	 * Exponential function.
	 *
	 * Computes exp(x), function result is nearly rounded. It will be correctly
	 * rounded to the theoretical value for 99.9% of input values, otherwise it
	 * will have a 1 UPL error.
	 *
	 * Method: Lookup intVal = exp(int(x)) Lookup fracVal = exp(int(x-int(x) /
	 * 1024.0) * 1024.0 ); Compute z as the exponential of the remaining bits by
	 * a polynomial minus one exp(x) = intVal * fracVal * (1 + z)
	 *
	 * Accuracy: Calculation is done with 63 bits of precision, so result should
	 * be correctly rounded for 99.9% of input values, with less than 1 ULP
	 * error otherwise.
	 *
	 * @param x
	 *            a double
	 * @return double e<sup>x</sup>
	 */
	public static double exp(double x) {
		return exp(x, 0.0, null);
	}

	/**
	 * Internal helper method for exponential function.
	 * 
	 * @param x
	 *            original argument of the exponential function
	 * @param extra
	 *            extra bits of precision on input (To Be Confirmed)
	 * @param hiPrec
	 *            extra bits of precision on output (To Be Confirmed)
	 * @return exp(x)
	 */
	private static double exp(double x, double extra, double[] hiPrec) {
		double intPartA;
		double intPartB;
		int intVal;

		/*
		 * Lookup exp(floor(x)). intPartA will have the upper 22 bits, intPartB
		 * will have the lower 52 bits.
		 */
		if (x < 0.0) {
			intVal = (int) -x;

			if (intVal > 746) {
				if (hiPrec != null) {
					hiPrec[0] = 0.0;
					hiPrec[1] = 0.0;
				}
				return 0.0;
			}

			if (intVal > 709) {
				/* This will produce a subnormal output */
				final double result = exp(x + 40.19140625, extra, hiPrec) / 285040095144011776.0;
				if (hiPrec != null) {
					hiPrec[0] /= 285040095144011776.0;
					hiPrec[1] /= 285040095144011776.0;
				}
				return result;
			}

			if (intVal == 709) {
				/* exp(1.494140625) is nearly a machine number... */
				final double result = exp(x + 1.494140625, extra, hiPrec) / 4.455505956692756620;
				if (hiPrec != null) {
					hiPrec[0] /= 4.455505956692756620;
					hiPrec[1] /= 4.455505956692756620;
				}
				return result;
			}

			intVal++;

			intPartA = EXP_INT_TABLE_A[750 - intVal];
			intPartB = EXP_INT_TABLE_B[750 - intVal];

			intVal = -intVal;
		} else {
			intVal = (int) x;

			if (intVal > 709) {
				if (hiPrec != null) {
					hiPrec[0] = Double.POSITIVE_INFINITY;
					hiPrec[1] = 0.0;
				}
				return Double.POSITIVE_INFINITY;
			}

			intPartA = EXP_INT_TABLE_A[750 + intVal];
			intPartB = EXP_INT_TABLE_B[750 + intVal];
		}

		/*
		 * Get the fractional part of x, find the greatest multiple of 2^-10
		 * less than x and look up the exp function of it. fracPartA will have
		 * the upper 22 bits, fracPartB the lower 52 bits.
		 */
		final int intFrac = (int) ((x - intVal) * 1024.0);
		final double fracPartA = EXP_FRAC_TABLE_A[intFrac];
		final double fracPartB = EXP_FRAC_TABLE_B[intFrac];

		/*
		 * epsilon is the difference in x from the nearest multiple of 2^-10. It
		 * has a value in the range 0 <= epsilon < 2^-10. Do the subtraction
		 * from x as the last step to avoid possible loss of percison.
		 */
		final double epsilon = x - (intVal + intFrac / 1024.0);

		/*
		 * Compute z = exp(epsilon) - 1.0 via a minimax polynomial. z has full
		 * double precision (52 bits). Since z < 2^-10, we will have 62 bits of
		 * precision when combined with the contant 1. This will be used in the
		 * last addition below to get proper rounding.
		 */

		/*
		 * Remez generated polynomial. Converges on the interval [0, 2^-10],
		 * error is less than 0.5 ULP
		 */
		double z = 0.04168701738764507;
		z = z * epsilon + 0.1666666505023083;
		z = z * epsilon + 0.5000000000042687;
		z = z * epsilon + 1.0;
		z = z * epsilon + -3.940510424527919E-20;

		/*
		 * Compute (intPartA+intPartB) * (fracPartA+fracPartB) by binomial
		 * expansion. tempA is exact since intPartA and intPartB only have 22
		 * bits each. tempB will have 52 bits of precision.
		 */
		double tempA = intPartA * fracPartA;
		double tempB = intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB;

		/*
		 * Compute the result. (1+z)(tempA+tempB). Order of operations is
		 * important. For accuracy add by increasing size. tempA is exact and
		 * much larger than the others. If there are extra bits specified from
		 * the pow() function, use them.
		 */
		final double tempC = tempB + tempA;
		final double result;
		if (extra != 0.0) {
			result = tempC * extra * z + tempC * extra + tempC * z + tempB + tempA;
		} else {
			result = tempC * z + tempB + tempA;
		}

		if (hiPrec != null) {
			// If requesting high precision
			hiPrec[0] = tempA;
			hiPrec[1] = tempC * extra * z + tempC * extra + tempC * z + tempB;
		}

		return result;
	}

	/**
	 * Compute exp(x) - 1
	 * 
	 * @param x
	 *            number to compute shifted exponential
	 * @return exp(x) - 1
	 */
	public static double expm1(double x) {
		return expm1(x, null);
	}

	/**
	 * Internal helper method for expm1
	 * 
	 * @param x
	 *            number to compute shifted exponential
	 * @param hiPrecOut
	 *            receive high precision result for -1.0 < x < 1.0
	 * @return exp(x) - 1
	 */
	private static double expm1(double x, double hiPrecOut[]) {
		if (x != x || x == 0.0) { // NaN or zero
			return x;
		}

		if (x <= -1.0 || x >= 1.0) {
			// If not between +/- 1.0
			// return exp(x) - 1.0;
			double hiPrec[] = new double[2];
			exp(x, 0.0, hiPrec);
			if (x > 0.0) {
				return -1.0 + hiPrec[0] + hiPrec[1];
			} else {
				final double ra = -1.0 + hiPrec[0];
				double rb = -(ra + 1.0 - hiPrec[0]);
				rb += hiPrec[1];
				return ra + rb;
			}
		}

		double baseA;
		double baseB;
		double epsilon;
		boolean negative = false;

		if (x < 0.0) {
			x = -x;
			negative = true;
		}

		{
			int intFrac = (int) (x * 1024.0);
			double tempA = EXP_FRAC_TABLE_A[intFrac] - 1.0;
			double tempB = EXP_FRAC_TABLE_B[intFrac];

			double temp = tempA + tempB;
			tempB = -(temp - tempA - tempB);
			tempA = temp;

			temp = tempA * HEX_40000000;
			baseA = tempA + temp - temp;
			baseB = tempB + (tempA - baseA);

			epsilon = x - intFrac / 1024.0;
		}

		/* Compute expm1(epsilon) */
		double zb = 0.008336750013465571;
		zb = zb * epsilon + 0.041666663879186654;
		zb = zb * epsilon + 0.16666666666745392;
		zb = zb * epsilon + 0.49999999999999994;
		zb = zb * epsilon;
		zb = zb * epsilon;

		double za = epsilon;
		double temp = za + zb;
		zb = -(temp - za - zb);
		za = temp;

		temp = za * HEX_40000000;
		temp = za + temp - temp;
		zb += za - temp;
		za = temp;

		/*
		 * Combine the parts. expm1(a+b) = expm1(a) + expm1(b) +
		 * expm1(a)*expm1(b)
		 */
		double ya = za * baseA;
		// double yb = za*baseB + zb*baseA + zb*baseB;
		temp = ya + za * baseB;
		double yb = -(temp - ya - za * baseB);
		ya = temp;

		temp = ya + zb * baseA;
		yb += -(temp - ya - zb * baseA);
		ya = temp;

		temp = ya + zb * baseB;
		yb += -(temp - ya - zb * baseB);
		ya = temp;

		// ya = ya + za + baseA;
		// yb = yb + zb + baseB;
		temp = ya + baseA;
		yb += -(temp - baseA - ya);
		ya = temp;

		temp = ya + za;
		// yb += (ya > za) ? -(temp - ya - za) : -(temp - za - ya);
		yb += -(temp - ya - za);
		ya = temp;

		temp = ya + baseB;
		// yb += (ya > baseB) ? -(temp - ya - baseB) : -(temp - baseB - ya);
		yb += -(temp - ya - baseB);
		ya = temp;

		temp = ya + zb;
		// yb += (ya > zb) ? -(temp - ya - zb) : -(temp - zb - ya);
		yb += -(temp - ya - zb);
		ya = temp;

		if (negative) {
			/* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
			double denom = 1.0 + ya;
			double denomr = 1.0 / denom;
			double denomb = -(denom - 1.0 - ya) + yb;
			double ratio = ya * denomr;
			temp = ratio * HEX_40000000;
			final double ra = ratio + temp - temp;
			double rb = ratio - ra;

			temp = denom * HEX_40000000;
			za = denom + temp - temp;
			zb = denom - za;

			rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;

			// f(x) = x/1+x
			// Compute f'(x)
			// Product rule: d(uv) = du*v + u*dv
			// Chain rule: d(f(g(x)) = f'(g(x))*f(g'(x))
			// d(1/x) = -1/(x*x)
			// d(1/1+x) = -1/( (1+x)^2) * 1 = -1/((1+x)*(1+x))
			// d(x/1+x) = -x/((1+x)(1+x)) + 1/1+x = 1 / ((1+x)(1+x))

			// Adjust for yb
			rb += yb * denomr; // numerator
			rb += -ya * denomb * denomr * denomr; // denominator

			// negate
			ya = -ra;
			yb = -rb;
		}

		if (hiPrecOut != null) {
			hiPrecOut[0] = ya;
			hiPrecOut[1] = yb;
		}

		return ya + yb;
	}

	/**
	 * For x between 0 and 1, returns exp(x), uses extended precision
	 * 
	 * @param x
	 *            argument of exponential
	 * @param result
	 *            placeholder where to place exp(x) split in two terms for extra
	 *            precision (i.e. exp(x) = result[0] + result[1]
	 * @return exp(x)
	 */
	private static double slowexp(final double x, final double result[]) {
		final double xs[] = new double[2];
		final double ys[] = new double[2];
		final double facts[] = new double[2];
		final double as[] = new double[2];
		split(x, xs);
		ys[0] = ys[1] = 0.0;

		for (int i = 19; i >= 0; i--) {
			splitMult(xs, ys, as);
			ys[0] = as[0];
			ys[1] = as[1];

			split(FACT[i], as);
			splitReciprocal(as, facts);

			splitAdd(ys, facts, as);
			ys[0] = as[0];
			ys[1] = as[1];
		}

		if (result != null) {
			result[0] = ys[0];
			result[1] = ys[1];
		}

		return ys[0] + ys[1];
	}

	/**
	 * Compute split[0], split[1] such that their sum is equal to d, and
	 * split[0] has its 30 least significant bits as zero.
	 * 
	 * @param d
	 *            number to split
	 * @param split
	 *            placeholder where to place the result
	 */
	private static void split(final double d, final double split[]) {
		if (d < 8e298 && d > -8e298) {
			final double a = d * HEX_40000000;
			split[0] = (d + a) - a;
			split[1] = d - split[0];
		} else {
			final double a = d * 9.31322574615478515625E-10;
			split[0] = (d + a - d) * HEX_40000000;
			split[1] = d - split[0];
		}
	}

	/**
	 * Recompute a split.
	 * 
	 * @param a
	 *            input/out array containing the split, changed on output
	 */
	private static void resplit(final double a[]) {
		final double c = a[0] + a[1];
		final double d = -(c - a[0] - a[1]);

		if (c < 8e298 && c > -8e298) {
			double z = c * HEX_40000000;
			a[0] = (c + z) - z;
			a[1] = c - a[0] + d;
		} else {
			double z = c * 9.31322574615478515625E-10;
			a[0] = (c + z - c) * HEX_40000000;
			a[1] = c - a[0] + d;
		}
	}

	/**
	 * Multiply two numbers in split form.
	 * 
	 * @param a
	 *            first term of multiplication
	 * @param b
	 *            second term of multiplication
	 * @param ans
	 *            placeholder where to put the result
	 */
	private static void splitMult(double a[], double b[], double ans[]) {
		ans[0] = a[0] * b[0];
		ans[1] = a[0] * b[1] + a[1] * b[0] + a[1] * b[1];

		/* Resplit */
		resplit(ans);
	}

	/**
	 * Add two numbers in split form.
	 * 
	 * @param a
	 *            first term of addition
	 * @param b
	 *            second term of addition
	 * @param ans
	 *            placeholder where to put the result
	 */
	private static void splitAdd(final double a[], final double b[], final double ans[]) {
		ans[0] = a[0] + b[0];
		ans[1] = a[1] + b[1];

		resplit(ans);
	}

	/**
	 * Compute the reciprocal of in. Use the following algorithm. in = c + d.
	 * want to find x + y such that x+y = 1/(c+d) and x is much larger than y
	 * and x has several zero bits on the right.
	 *
	 * Set b = 1/(2^22), a = 1 - b. Thus (a+b) = 1. Use following identity to
	 * compute (a+b)/(c+d)
	 *
	 * (a+b)/(c+d) = a/c + (bc - ad) / (c^2 + cd) set x = a/c and y = (bc - ad)
	 * / (c^2 + cd) This will be close to the right answer, but there will be
	 * some rounding in the calculation of X. So by carefully computing 1 -
	 * (c+d)(x+y) we can compute an error and add that back in. This is done
	 * carefully so that terms of similar size are subtracted first.
	 * 
	 * @param in
	 *            initial number, in split form
	 * @param result
	 *            placeholder where to put the result
	 */
	private static void splitReciprocal(final double in[], final double result[]) {
		final double b = 1.0 / 4194304.0;
		final double a = 1.0 - b;

		if (in[0] == 0.0) {
			in[0] = in[1];
			in[1] = 0.0;
		}

		result[0] = a / in[0];
		result[1] = (b * in[0] - a * in[1]) / (in[0] * in[0] + in[0] * in[1]);

		if (result[1] != result[1]) { // can happen if result[1] is NAN
			result[1] = 0.0;
		}

		/* Resplit */
		resplit(result);

		for (int i = 0; i < 2; i++) {
			/* this may be overkill, probably once is enough */
			double err = 1.0 - result[0] * in[0] - result[0] * in[1] - result[1] * in[0] - result[1] * in[1];
			/* err = 1.0 - err; */
			err = err * (result[0] + result[1]);
			/* printf("err = %16e\n", err); */
			result[1] += err;
		}
	}

	/**
	 * Compute (a[0] + a[1]) * (b[0] + b[1]) in extended precision.
	 * 
	 * @param a
	 *            first term of the multiplication
	 * @param b
	 *            second term of the multiplication
	 * @param result
	 *            placeholder where to put the result
	 */
	private static void quadMult(final double a[], final double b[], final double result[]) {
		final double xs[] = new double[2];
		final double ys[] = new double[2];
		final double zs[] = new double[2];

		/* a[0] * b[0] */
		split(a[0], xs);
		split(b[0], ys);
		splitMult(xs, ys, zs);

		result[0] = zs[0];
		result[1] = zs[1];

		/* a[0] * b[1] */
		split(b[1], ys);
		splitMult(xs, ys, zs);

		double tmp = result[0] + zs[0];
		result[1] = result[1] - (tmp - result[0] - zs[0]);
		result[0] = tmp;
		tmp = result[0] + zs[1];
		result[1] = result[1] - (tmp - result[0] - zs[1]);
		result[0] = tmp;

		/* a[1] * b[0] */
		split(a[1], xs);
		split(b[0], ys);
		splitMult(xs, ys, zs);

		tmp = result[0] + zs[0];
		result[1] = result[1] - (tmp - result[0] - zs[0]);
		result[0] = tmp;
		tmp = result[0] + zs[1];
		result[1] = result[1] - (tmp - result[0] - zs[1]);
		result[0] = tmp;

		/* a[1] * b[0] */
		split(a[1], xs);
		split(b[1], ys);
		splitMult(xs, ys, zs);

		tmp = result[0] + zs[0];
		result[1] = result[1] - (tmp - result[0] - zs[0]);
		result[0] = tmp;
		tmp = result[0] + zs[1];
		result[1] = result[1] - (tmp - result[0] - zs[1]);
		result[0] = tmp;
	}

	/**
	 * Compute exp(p) for a integer p in extended precision.
	 * 
	 * @param p
	 *            integer whose exponential is requested
	 * @param result
	 *            placeholder where to put the result in extended precision
	 * @return exp(p) in standard precision (equal to result[0] + result[1])
	 */
	private static double expint(int p, final double result[]) {
		// double x = M_E;
		final double xs[] = new double[2];
		final double as[] = new double[2];
		final double ys[] = new double[2];
		// split(x, xs);
		// xs[1] = (double)(2.7182818284590452353602874713526625L - xs[0]);
		// xs[0] = 2.71827697753906250000;
		// xs[1] = 4.85091998273542816811e-06;
		// xs[0] = Double.longBitsToDouble(0x4005bf0800000000L);
		// xs[1] = Double.longBitsToDouble(0x3ed458a2bb4a9b00L);

		/* E */
		xs[0] = 2.718281828459045;
		xs[1] = 1.4456468917292502E-16;

		split(1.0, ys);

		while (p > 0) {
			if ((p & 1) != 0) {
				quadMult(ys, xs, as);
				ys[0] = as[0];
				ys[1] = as[1];
			}

			quadMult(xs, xs, as);
			xs[0] = as[0];
			xs[1] = as[1];

			p >>= 1;
		}

		if (result != null) {
			result[0] = ys[0];
			result[1] = ys[1];

			resplit(result);
		}

		return ys[0] + ys[1];
	}

	/**
	 * Natural logarithm.
	 *
	 * @param x
	 *            a double
	 * @return log(x)
	 */
	public static double log(final double x) {
		return log(x, null);
	}

	/**
	 * Internal helper method for natural logarithm function.
	 * 
	 * @param x
	 *            original argument of the natural logarithm function
	 * @param hiPrec
	 *            extra bits of precision on output (To Be Confirmed)
	 * @return log(x)
	 */
	private static double log(final double x, final double[] hiPrec) {
		if (x == 0) { // Handle special case of +0/-0
			return Double.NEGATIVE_INFINITY;
		}
		long bits = Double.doubleToLongBits(x);

		/* Handle special cases of negative input, and NaN */
		if ((bits & 0x8000000000000000L) != 0 || x != x) {
			if (x != 0.0) {
				if (hiPrec != null) {
					hiPrec[0] = Double.NaN;
				}

				return Double.NaN;
			}
		}

		/* Handle special cases of Positive infinity. */
		if (x == Double.POSITIVE_INFINITY) {
			if (hiPrec != null) {
				hiPrec[0] = Double.POSITIVE_INFINITY;
			}

			return Double.POSITIVE_INFINITY;
		}

		/* Extract the exponent */
		int exp = (int) (bits >> 52) - 1023;

		if ((bits & 0x7ff0000000000000L) == 0) {
			// Subnormal!
			if (x == 0) {
				// Zero
				if (hiPrec != null) {
					hiPrec[0] = Double.NEGATIVE_INFINITY;
				}

				return Double.NEGATIVE_INFINITY;
			}

			/* Normalize the subnormal number. */
			bits <<= 1;
			while ((bits & 0x0010000000000000L) == 0) {
				exp--;
				bits <<= 1;
			}
		}

		if (exp == -1 || exp == 0) {
			if (x < 1.01 && x > 0.99 && hiPrec == null) {
				/*
				 * The normal method doesn't work well in the range [0.99,
				 * 1.01], so call do a straight polynomial expansion in higer
				 * precision.
				 */

				/* Compute x - 1.0 and split it */
				double xa = x - 1.0;
				double xb = xa - x + 1.0;
				double tmp = xa * HEX_40000000;
				double aa = xa + tmp - tmp;
				double ab = xa - aa;
				xa = aa;
				xb = ab;

				double ya = LN_QUICK_COEF[LN_QUICK_COEF.length - 1][0];
				double yb = LN_QUICK_COEF[LN_QUICK_COEF.length - 1][1];

				for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) {
					/* Multiply a = y * x */
					aa = ya * xa;
					ab = ya * xb + yb * xa + yb * xb;
					/* split, so now y = a */
					tmp = aa * HEX_40000000;
					ya = aa + tmp - tmp;
					yb = aa - ya + ab;

					/* Add a = y + lnQuickCoef */
					aa = ya + LN_QUICK_COEF[i][0];
					ab = yb + LN_QUICK_COEF[i][1];
					/* Split y = a */
					tmp = aa * HEX_40000000;
					ya = aa + tmp - tmp;
					yb = aa - ya + ab;
				}

				/* Multiply a = y * x */
				aa = ya * xa;
				ab = ya * xb + yb * xa + yb * xb;
				/* split, so now y = a */
				tmp = aa * HEX_40000000;
				ya = aa + tmp - tmp;
				yb = aa - ya + ab;

				return ya + yb;
			}
		}

		// lnm is a log of a number in the range of 1.0 - 2.0, so 0 <= lnm <
		// ln(2)
		double lnm[] = LN_MANT[(int) ((bits & 0x000ffc0000000000L) >> 42)];

		/*
		 * double epsilon = x / Double.longBitsToDouble(bits &
		 * 0xfffffc0000000000L);
		 * 
		 * epsilon -= 1.0;
		 */

		// y is the most significant 10 bits of the mantissa
		// double y = Double.longBitsToDouble(bits & 0xfffffc0000000000L);
		// double epsilon = (x - y) / y;
		double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));

		double lnza = 0.0;
		double lnzb = 0.0;

		if (hiPrec != null) {
			/* split epsilon -> x */
			double tmp = epsilon * HEX_40000000;
			double aa = epsilon + tmp - tmp;
			double ab = epsilon - aa;
			double xa = aa;
			double xb = ab;

			/* Need a more accurate epsilon, so adjust the division. */
			double numer = bits & 0x3ffffffffffL;
			double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
			aa = numer - xa * denom - xb * denom;
			xb += aa / denom;

			/* Remez polynomial evaluation */
			double ya = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length - 1][0];
			double yb = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length - 1][1];

			for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) {
				/* Multiply a = y * x */
				aa = ya * xa;
				ab = ya * xb + yb * xa + yb * xb;
				/* split, so now y = a */
				tmp = aa * HEX_40000000;
				ya = aa + tmp - tmp;
				yb = aa - ya + ab;

				/* Add a = y + lnHiPrecCoef */
				aa = ya + LN_HI_PREC_COEF[i][0];
				ab = yb + LN_HI_PREC_COEF[i][1];
				/* Split y = a */
				tmp = aa * HEX_40000000;
				ya = aa + tmp - tmp;
				yb = aa - ya + ab;
			}

			/* Multiply a = y * x */
			aa = ya * xa;
			ab = ya * xb + yb * xa + yb * xb;

			/* split, so now lnz = a */
			/*
			 * tmp = aa * 1073741824.0; lnza = aa + tmp - tmp; lnzb = aa - lnza
			 * + ab;
			 */
			lnza = aa + ab;
			lnzb = -(lnza - aa - ab);
		} else {
			/*
			 * High precision not required. Eval Remez polynomial using standard
			 * double precision
			 */
			lnza = -0.16624882440418567;
			lnza = lnza * epsilon + 0.19999954120254515;
			lnza = lnza * epsilon + -0.2499999997677497;
			lnza = lnza * epsilon + 0.3333333333332802;
			lnza = lnza * epsilon + -0.5;
			lnza = lnza * epsilon + 1.0;
			lnza = lnza * epsilon;
		}

		/*
		 * Relative sizes: lnzb [0, 2.33E-10] lnm[1] [0, 1.17E-7] ln2B*exp [0,
		 * 1.12E-4] lnza [0, 9.7E-4] lnm[0] [0, 0.692] ln2A*exp [0, 709]
		 */

		/*
		 * Compute the following sum: lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] +
		 * ln2A*exp;
		 */

		// return lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
		double a = LN_2_A * exp;
		double b = 0.0;
		double c = a + lnm[0];
		double d = -(c - a - lnm[0]);
		a = c;
		b = b + d;

		c = a + lnza;
		d = -(c - a - lnza);
		a = c;
		b = b + d;

		c = a + LN_2_B * exp;
		d = -(c - a - LN_2_B * exp);
		a = c;
		b = b + d;

		c = a + lnm[1];
		d = -(c - a - lnm[1]);
		a = c;
		b = b + d;

		c = a + lnzb;
		d = -(c - a - lnzb);
		a = c;
		b = b + d;

		if (hiPrec != null) {
			hiPrec[0] = a;
			hiPrec[1] = b;
		}

		return a + b;
	}

	/**
	 * Compute log(1 + x).
	 * 
	 * @param x
	 *            a number
	 * @return log(1 + x)
	 */
	public static double log1p(final double x) {
		double xpa = 1.0 + x;
		double xpb = -(xpa - 1.0 - x);

		if (x == -1) {
			return x / 0.0; // -Infinity
		}

		if (x > 0 && 1 / x == 0) { // x = Infinity
			return x;
		}

		if (x > 1e-6 || x < -1e-6) {
			double hiPrec[] = new double[2];

			final double lores = log(xpa, hiPrec);
			if (Double.isInfinite(lores)) { // don't allow this to be converted
											// to NaN
				return lores;
			}

			/* Do a taylor series expansion around xpa */
			/* f(x+y) = f(x) + f'(x)*y + f''(x)/2 y^2 */
			double fx1 = xpb / xpa;

			double epsilon = 0.5 * fx1 + 1.0;
			epsilon = epsilon * fx1;

			return epsilon + hiPrec[1] + hiPrec[0];
		}

		/* Value is small |x| < 1e6, do a Taylor series centered on 1.0 */
		double y = x * 0.333333333333333 - 0.5;
		y = y * x + 1.0;
		y = y * x;

		return y;
	}

	/**
	 * Compute the base 10 logarithm.
	 * 
	 * @param x
	 *            a number
	 * @return log10(x)
	 */
	public static double log10(final double x) {
		final double hiPrec[] = new double[2];

		final double lores = log(x, hiPrec);
		if (Double.isInfinite(lores)) { // don't allow this to be converted to
										// NaN
			return lores;
		}

		final double tmp = hiPrec[0] * HEX_40000000;
		final double lna = hiPrec[0] + tmp - tmp;
		final double lnb = hiPrec[0] - lna + hiPrec[1];

		final double rln10a = 0.4342944622039795;
		final double rln10b = 1.9699272335463627E-8;

		return rln10b * lnb + rln10b * lna + rln10a * lnb + rln10a * lna;
	}

	/**
	 * Power function. Compute x^y.
	 *
	 * @param x
	 *            a double
	 * @param y
	 *            a double
	 * @return double
	 */
	public static double pow(double x, double y) {
		final double lns[] = new double[2];

		if (y == 0.0) {
			return 1.0;
		}

		if (x != x) { // X is NaN
			return x;
		}

		if (x == 0) {
			long bits = Double.doubleToLongBits(x);
			if ((bits & 0x8000000000000000L) != 0) {
				// -zero
				long yi = (long) y;

				if (y < 0 && y == yi && (yi & 1) == 1) {
					return Double.NEGATIVE_INFINITY;
				}

				if (y < 0 && y == yi && (yi & 1) == 1) {
					return -0.0;
				}

				if (y > 0 && y == yi && (yi & 1) == 1) {
					return -0.0;
				}
			}

			if (y < 0) {
				return Double.POSITIVE_INFINITY;
			}
			if (y > 0) {
				return 0.0;
			}

			return Double.NaN;
		}

		if (x == Double.POSITIVE_INFINITY) {
			if (y != y) { // y is NaN
				return y;
			}
			if (y < 0.0) {
				return 0.0;
			} else {
				return Double.POSITIVE_INFINITY;
			}
		}

		if (y == Double.POSITIVE_INFINITY) {
			if (x * x == 1.0)
				return Double.NaN;

			if (x * x > 1.0) {
				return Double.POSITIVE_INFINITY;
			} else {
				return 0.0;
			}
		}

		if (x == Double.NEGATIVE_INFINITY) {
			if (y != y) { // y is NaN
				return y;
			}

			if (y < 0) {
				long yi = (long) y;
				if (y == yi && (yi & 1) == 1) {
					return -0.0;
				}

				return 0.0;
			}

			if (y > 0) {
				long yi = (long) y;
				if (y == yi && (yi & 1) == 1) {
					return Double.NEGATIVE_INFINITY;
				}

				return Double.POSITIVE_INFINITY;
			}
		}

		if (y == Double.NEGATIVE_INFINITY) {

			if (x * x == 1.0) {
				return Double.NaN;
			}

			if (x * x < 1.0) {
				return Double.POSITIVE_INFINITY;
			} else {
				return 0.0;
			}
		}

		/* Handle special case x<0 */
		if (x < 0) {
			// y is an even integer in this case
			if (y >= TWO_POWER_52 || y <= -TWO_POWER_52) {
				return pow(-x, y);
			}

			if (y == (long) y) {
				// If y is an integer
				return ((long) y & 1) == 0 ? pow(-x, y) : -pow(-x, y);
			} else {
				return Double.NaN;
			}
		}

		/* Split y into ya and yb such that y = ya+yb */
		double ya;
		double yb;
		if (y < 8e298 && y > -8e298) {
			double tmp1 = y * HEX_40000000;
			ya = y + tmp1 - tmp1;
			yb = y - ya;
		} else {
			double tmp1 = y * 9.31322574615478515625E-10;
			double tmp2 = tmp1 * 9.31322574615478515625E-10;
			ya = (tmp1 + tmp2 - tmp1) * HEX_40000000 * HEX_40000000;
			yb = y - ya;
		}

		/* Compute ln(x) */
		final double lores = log(x, lns);
		if (Double.isInfinite(lores)) { // don't allow this to be converted to
										// NaN
			return lores;
		}

		double lna = lns[0];
		double lnb = lns[1];

		/* resplit lns */
		double tmp1 = lna * HEX_40000000;
		double tmp2 = lna + tmp1 - tmp1;
		lnb += lna - tmp2;
		lna = tmp2;

		// y*ln(x) = (aa+ab)
		final double aa = lna * ya;
		final double ab = lna * yb + lnb * ya + lnb * yb;

		lna = aa + ab;
		lnb = -(lna - aa - ab);

		double z = 1.0 / 120.0;
		z = z * lnb + (1.0 / 24.0);
		z = z * lnb + (1.0 / 6.0);
		z = z * lnb + 0.5;
		z = z * lnb + 1.0;
		z = z * lnb;

		final double result = exp(lna, z, null);
		// result = result + result * z;
		return result;
	}

	/**
	 * xi in the range of [1, 2]. 3 5 7 x+1 / x x x \ ln ----- = 2 * | x + ----
	 * + ---- + ---- + ... | 1-x \ 3 5 7 /
	 *
	 * So, compute a Remez approximation of the following function
	 *
	 * ln ((sqrt(x)+1)/(1-sqrt(x))) / x
	 *
	 * This will be an even function with only positive coefficents. x is in the
	 * range [0 - 1/3].
	 *
	 * Transform xi for input to the above function by setting x =
	 * (xi-1)/(xi+1). Input to the polynomial is x^2, then the result is
	 * multiplied by x.
	 * 
	 * @param xi
	 *            number from which log is requested
	 * @return log(xi)
	 */
	private static double[] slowLog(double xi) {
		double x[] = new double[2];
		double x2[] = new double[2];
		double y[] = new double[2];
		double a[] = new double[2];

		split(xi, x);

		/* Set X = (x-1)/(x+1) */
		x[0] += 1.0;
		resplit(x);
		splitReciprocal(x, a);
		x[0] -= 2.0;
		resplit(x);
		splitMult(x, a, y);
		x[0] = y[0];
		x[1] = y[1];

		/* Square X -> X2 */
		splitMult(x, x, x2);

		// x[0] -= 1.0;
		// resplit(x);

		y[0] = LN_SPLIT_COEF[LN_SPLIT_COEF.length - 1][0];
		y[1] = LN_SPLIT_COEF[LN_SPLIT_COEF.length - 1][1];

		for (int i = LN_SPLIT_COEF.length - 2; i >= 0; i--) {
			splitMult(y, x2, a);
			y[0] = a[0];
			y[1] = a[1];
			splitAdd(y, LN_SPLIT_COEF[i], a);
			y[0] = a[0];
			y[1] = a[1];
		}

		splitMult(y, x, a);
		y[0] = a[0];
		y[1] = a[1];

		return y;
	}

	/**
	 * For x between 0 and pi/4 compute sine.
	 * 
	 * @param x
	 *            number from which sine is requested
	 * @param result
	 *            placeholder where to put the result in extended precision
	 * @return sin(x)
	 */
	private static double slowSin(final double x, final double result[]) {
		final double xs[] = new double[2];
		final double ys[] = new double[2];
		final double facts[] = new double[2];
		final double as[] = new double[2];
		split(x, xs);
		ys[0] = ys[1] = 0.0;

		for (int i = 19; i >= 0; i--) {
			splitMult(xs, ys, as);
			ys[0] = as[0];
			ys[1] = as[1];

			if ((i & 1) == 0) {
				continue;
			}

			split(FACT[i], as);
			splitReciprocal(as, facts);

			if ((i & 2) != 0) {
				facts[0] = -facts[0];
				facts[1] = -facts[1];
			}

			splitAdd(ys, facts, as);
			ys[0] = as[0];
			ys[1] = as[1];
		}

		if (result != null) {
			result[0] = ys[0];
			result[1] = ys[1];
		}

		return ys[0] + ys[1];
	}

	/**
	 * For x between 0 and pi/4 compute cosine
	 * 
	 * @param x
	 *            number from which cosine is requested
	 * @param result
	 *            placeholder where to put the result in extended precision
	 * @return cos(x)
	 */
	private static double slowCos(final double x, final double result[]) {

		final double xs[] = new double[2];
		final double ys[] = new double[2];
		final double facts[] = new double[2];
		final double as[] = new double[2];
		split(x, xs);
		ys[0] = ys[1] = 0.0;

		for (int i = 19; i >= 0; i--) {
			splitMult(xs, ys, as);
			ys[0] = as[0];
			ys[1] = as[1];

			if ((i & 1) != 0) {
				continue;
			}

			split(FACT[i], as);
			splitReciprocal(as, facts);

			if ((i & 2) != 0) {
				facts[0] = -facts[0];
				facts[1] = -facts[1];
			}

			splitAdd(ys, facts, as);
			ys[0] = as[0];
			ys[1] = as[1];
		}

		if (result != null) {
			result[0] = ys[0];
			result[1] = ys[1];
		}

		return ys[0] + ys[1];
	}

	/**
	 * Build the sine and cosine tables.
	 */
	private static void buildSinCosTables() {
		final double result[] = new double[2];

		/* Use taylor series for 0 <= x <= 6/8 */
		for (int i = 0; i < 7; i++) {
			double x = i / 8.0;

			slowSin(x, result);
			SINE_TABLE_A[i] = result[0];
			SINE_TABLE_B[i] = result[1];

			slowCos(x, result);
			COSINE_TABLE_A[i] = result[0];
			COSINE_TABLE_B[i] = result[1];
		}

		/*
		 * Use angle addition formula to complete table to 13/8, just beyond
		 * pi/2
		 */
		for (int i = 7; i < 14; i++) {
			double xs[] = new double[2];
			double ys[] = new double[2];
			double as[] = new double[2];
			double bs[] = new double[2];
			double temps[] = new double[2];

			if ((i & 1) == 0) {
				// Even, use double angle
				xs[0] = SINE_TABLE_A[i / 2];
				xs[1] = SINE_TABLE_B[i / 2];
				ys[0] = COSINE_TABLE_A[i / 2];
				ys[1] = COSINE_TABLE_B[i / 2];

				/* compute sine */
				splitMult(xs, ys, result);
				SINE_TABLE_A[i] = result[0] * 2.0;
				SINE_TABLE_B[i] = result[1] * 2.0;

				/* Compute cosine */
				splitMult(ys, ys, as);
				splitMult(xs, xs, temps);
				temps[0] = -temps[0];
				temps[1] = -temps[1];
				splitAdd(as, temps, result);
				COSINE_TABLE_A[i] = result[0];
				COSINE_TABLE_B[i] = result[1];
			} else {
				xs[0] = SINE_TABLE_A[i / 2];
				xs[1] = SINE_TABLE_B[i / 2];
				ys[0] = COSINE_TABLE_A[i / 2];
				ys[1] = COSINE_TABLE_B[i / 2];
				as[0] = SINE_TABLE_A[i / 2 + 1];
				as[1] = SINE_TABLE_B[i / 2 + 1];
				bs[0] = COSINE_TABLE_A[i / 2 + 1];
				bs[1] = COSINE_TABLE_B[i / 2 + 1];

				/* compute sine */
				splitMult(xs, bs, temps);
				splitMult(ys, as, result);
				splitAdd(result, temps, result);
				SINE_TABLE_A[i] = result[0];
				SINE_TABLE_B[i] = result[1];

				/* Compute cosine */
				splitMult(ys, bs, result);
				splitMult(xs, as, temps);
				temps[0] = -temps[0];
				temps[1] = -temps[1];
				splitAdd(result, temps, result);
				COSINE_TABLE_A[i] = result[0];
				COSINE_TABLE_B[i] = result[1];
			}
		}

		/* Compute tangent = sine/cosine */
		for (int i = 0; i < 14; i++) {
			double xs[] = new double[2];
			double ys[] = new double[2];
			double as[] = new double[2];

			as[0] = COSINE_TABLE_A[i];
			as[1] = COSINE_TABLE_B[i];

			splitReciprocal(as, ys);

			xs[0] = SINE_TABLE_A[i];
			xs[1] = SINE_TABLE_B[i];

			splitMult(xs, ys, as);

			TANGENT_TABLE_A[i] = as[0];
			TANGENT_TABLE_B[i] = as[1];
		}

	}

	/**
	 * Computes sin(x) - x, where |x| < 1/16. Use a Remez polynomial
	 * approximation.
	 * 
	 * @param x
	 *            a number smaller than 1/16
	 * @return sin(x) - x
	 */
	private static double polySine(final double x) {
		double x2 = x * x;

		double p = 2.7553817452272217E-6;
		p = p * x2 + -1.9841269659586505E-4;
		p = p * x2 + 0.008333333333329196;
		p = p * x2 + -0.16666666666666666;
		// p *= x2;
		// p *= x;
		p = p * x2 * x;

		return p;
	}

	/**
	 * Computes cos(x) - 1, where |x| < 1/16. Use a Remez polynomial
	 * approximation.
	 * 
	 * @param x
	 *            a number smaller than 1/16
	 * @return cos(x) - 1
	 */
	private static double polyCosine(double x) {
		double x2 = x * x;

		double p = 2.479773539153719E-5;
		p = p * x2 + -0.0013888888689039883;
		p = p * x2 + 0.041666666666621166;
		p = p * x2 + -0.49999999999999994;
		p *= x2;

		return p;
	}

	/**
	 * Compute sine over the first quadrant (0 < x < pi/2). Use combination of
	 * table lookup and rational polynomial expansion.
	 * 
	 * @param xa
	 *            number from which sine is requested
	 * @param xb
	 *            extra bits for x (may be 0.0)
	 * @return sin(xa + xb)
	 */
	private static double sinQ(double xa, double xb) {
		int idx = (int) ((xa * 8.0) + 0.5);
		final double epsilon = xa - EIGHTHS[idx]; // idx*0.125;

		// Table lookups
		final double sintA = SINE_TABLE_A[idx];
		final double sintB = SINE_TABLE_B[idx];
		final double costA = COSINE_TABLE_A[idx];
		final double costB = COSINE_TABLE_B[idx];

		// Polynomial eval of sin(epsilon), cos(epsilon)
		double sinEpsA = epsilon;
		double sinEpsB = polySine(epsilon);
		final double cosEpsA = 1.0;
		final double cosEpsB = polyCosine(epsilon);

		// Split epsilon xa + xb = x
		final double temp = sinEpsA * HEX_40000000;
		double temp2 = (sinEpsA + temp) - temp;
		sinEpsB += sinEpsA - temp2;
		sinEpsA = temp2;

		/* Compute sin(x) by angle addition formula */
		double result;

		/*
		 * Compute the following sum:
		 * 
		 * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
		 * sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
		 * 
		 * Ranges of elements
		 * 
		 * xxxtA 0 PI/2 xxxtB -1.5e-9 1.5e-9 sinEpsA -0.0625 0.0625 sinEpsB
		 * -6e-11 6e-11 cosEpsA 1.0 cosEpsB 0 -0.0625
		 */

		// result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
		// sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;

		// result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
		// result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB *
		// sinEpsB;
		double a = 0;
		double b = 0;

		double t = sintA;
		double c = a + t;
		double d = -(c - a - t);
		a = c;
		b = b + d;

		t = costA * sinEpsA;
		c = a + t;
		d = -(c - a - t);
		a = c;
		b = b + d;

		b = b + sintA * cosEpsB + costA * sinEpsB;
		/*
		 * t = sintA*cosEpsB; c = a + t; d = -(c - a - t); a = c; b = b + d;
		 * 
		 * t = costA*sinEpsB; c = a + t; d = -(c - a - t); a = c; b = b + d;
		 */

		b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;
		/*
		 * t = sintB; c = a + t; d = -(c - a - t); a = c; b = b + d;
		 * 
		 * t = costB*sinEpsA; c = a + t; d = -(c - a - t); a = c; b = b + d;
		 * 
		 * t = sintB*cosEpsB; c = a + t; d = -(c - a - t); a = c; b = b + d;
		 * 
		 * t = costB*sinEpsB; c = a + t; d = -(c - a - t); a = c; b = b + d;
		 */

		if (xb != 0.0) {
			t = ((costA + costB) * (cosEpsA + cosEpsB) - (sintA + sintB) * (sinEpsA + sinEpsB)) * xb; // approximate
																										// cosine*xb
			c = a + t;
			d = -(c - a - t);
			a = c;
			b = b + d;
		}

		result = a + b;

		return result;
	}

	/**
	 * Compute cosine in the first quadrant by subtracting input from PI/2 and
	 * then calling sinQ. This is more accurate as the input approaches PI/2.
	 * 
	 * @param xa
	 *            number from which cosine is requested
	 * @param xb
	 *            extra bits for x (may be 0.0)
	 * @return cos(xa + xb)
	 */
	private static double cosQ(double xa, double xb) {
		final double pi2a = 1.5707963267948966;
		final double pi2b = 6.123233995736766E-17;

		final double a = pi2a - xa;
		double b = -(a - pi2a + xa);
		b += pi2b - xb;

		return sinQ(a, b);
	}

	/**
	 * Compute tangent (or cotangent) over the first quadrant. 0 < x < pi/2 Use
	 * combination of table lookup and rational polynomial expansion.
	 * 
	 * @param xa
	 *            number from which sine is requested
	 * @param xb
	 *            extra bits for x (may be 0.0)
	 * @param cotanFlag
	 *            if true, compute the cotangent instead of the tangent
	 * @return tan(xa+xb) (or cotangent, depending on cotanFlag)
	 */
	private static double tanQ(double xa, double xb, boolean cotanFlag) {

		int idx = (int) ((xa * 8.0) + 0.5);
		final double epsilon = xa - EIGHTHS[idx]; // idx*0.125;

		// Table lookups
		final double sintA = SINE_TABLE_A[idx];
		final double sintB = SINE_TABLE_B[idx];
		final double costA = COSINE_TABLE_A[idx];
		final double costB = COSINE_TABLE_B[idx];

		// Polynomial eval of sin(epsilon), cos(epsilon)
		double sinEpsA = epsilon;
		double sinEpsB = polySine(epsilon);
		final double cosEpsA = 1.0;
		final double cosEpsB = polyCosine(epsilon);

		// Split epsilon xa + xb = x
		double temp = sinEpsA * HEX_40000000;
		double temp2 = (sinEpsA + temp) - temp;
		sinEpsB += sinEpsA - temp2;
		sinEpsA = temp2;

		/* Compute sin(x) by angle addition formula */

		/*
		 * Compute the following sum:
		 * 
		 * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
		 * sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
		 * 
		 * Ranges of elements
		 * 
		 * xxxtA 0 PI/2 xxxtB -1.5e-9 1.5e-9 sinEpsA -0.0625 0.0625 sinEpsB
		 * -6e-11 6e-11 cosEpsA 1.0 cosEpsB 0 -0.0625
		 */

		// result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
		// sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;

		// result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
		// result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB *
		// sinEpsB;
		double a = 0;
		double b = 0;

		// Compute sine
		double t = sintA;
		double c = a + t;
		double d = -(c - a - t);
		a = c;
		b = b + d;

		t = costA * sinEpsA;
		c = a + t;
		d = -(c - a - t);
		a = c;
		b = b + d;

		b = b + sintA * cosEpsB + costA * sinEpsB;
		b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;

		double sina = a + b;
		double sinb = -(sina - a - b);

		// Compute cosine

		a = b = c = d = 0.0;

		t = costA * cosEpsA;
		c = a + t;
		d = -(c - a - t);
		a = c;
		b = b + d;

		t = -sintA * sinEpsA;
		c = a + t;
		d = -(c - a - t);
		a = c;
		b = b + d;

		b = b + costB * cosEpsA + costA * cosEpsB + costB * cosEpsB;
		b = b - (sintB * sinEpsA + sintA * sinEpsB + sintB * sinEpsB);

		double cosa = a + b;
		double cosb = -(cosa - a - b);

		if (cotanFlag) {
			double tmp;
			tmp = cosa;
			cosa = sina;
			sina = tmp;
			tmp = cosb;
			cosb = sinb;
			sinb = tmp;
		}

		/* estimate and correct, compute 1.0/(cosa+cosb) */
		/*
		 * double est = (sina+sinb)/(cosa+cosb); double err = (sina - cosa*est)
		 * + (sinb - cosb*est); est += err/(cosa+cosb); err = (sina - cosa*est)
		 * + (sinb - cosb*est);
		 */

		// f(x) = 1/x, f'(x) = -1/x^2

		double est = sina / cosa;

		/* Split the estimate to get more accurate read on division rounding */
		temp = est * HEX_40000000;
		double esta = (est + temp) - temp;
		double estb = est - esta;

		temp = cosa * HEX_40000000;
		double cosaa = (cosa + temp) - temp;
		double cosab = cosa - cosaa;

		// double err = (sina - est*cosa)/cosa; // Correction for division
		// rounding
		double err = (sina - esta * cosaa - esta * cosab - estb * cosaa - estb * cosab) / cosa; // Correction
																								// for
																								// division
																								// rounding
		err += sinb / cosa; // Change in est due to sinb
		err += -sina * cosb / cosa / cosa; // Change in est due to cosb

		if (xb != 0.0) {
			// tan' = 1 + tan^2 cot' = -(1 + cot^2)
			// Approximate impact of xb
			double xbadj = xb + est * est * xb;
			if (cotanFlag) {
				xbadj = -xbadj;
			}

			err += xbadj;
		}

		return est + err;
	}

	/**
	 * Reduce the input argument using the Payne and Hanek method. This is good
	 * for all inputs 0.0 < x < inf Output is remainder after dividing by PI/2
	 * The result array should contain 3 numbers. result[0] is the integer
	 * portion, so mod 4 this gives the quadrant. result[1] is the upper bits of
	 * the remainder result[2] is the lower bits of the remainder
	 *
	 * @param x
	 *            number to reduce
	 * @param result
	 *            placeholder where to put the result
	 */
	private static void reducePayneHanek(double x, double result[]) {
		/* Convert input double to bits */
		long inbits = Double.doubleToLongBits(x);
		int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;

		/* Convert to fixed point representation */
		inbits &= 0x000fffffffffffffL;
		inbits |= 0x0010000000000000L;

		/* Normalize input to be between 0.5 and 1.0 */
		exponent++;
		inbits <<= 11;

		/* Based on the exponent, get a shifted copy of recip2pi */
		long shpi0;
		long shpiA;
		long shpiB;
		int idx = exponent >> 6;
		int shift = exponent - (idx << 6);

		if (shift != 0) {
			shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx - 1] << shift);
			shpi0 |= RECIP_2PI[idx] >>> (64 - shift);
			shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx + 1] >>> (64 - shift));
			shpiB = (RECIP_2PI[idx + 1] << shift) | (RECIP_2PI[idx + 2] >>> (64 - shift));
		} else {
			shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx - 1];
			shpiA = RECIP_2PI[idx];
			shpiB = RECIP_2PI[idx + 1];
		}

		/* Multiply input by shpiA */
		long a = inbits >>> 32;
		long b = inbits & 0xffffffffL;

		long c = shpiA >>> 32;
		long d = shpiA & 0xffffffffL;

		long ac = a * c;
		long bd = b * d;
		long bc = b * c;
		long ad = a * d;

		long prodB = bd + (ad << 32);
		long prodA = ac + (ad >>> 32);

		boolean bita = (bd & 0x8000000000000000L) != 0;
		boolean bitb = (ad & 0x80000000L) != 0;
		boolean bitsum = (prodB & 0x8000000000000000L) != 0;

		/* Carry */
		if ((bita && bitb) || ((bita || bitb) && !bitsum)) {
			prodA++;
		}

		bita = (prodB & 0x8000000000000000L) != 0;
		bitb = (bc & 0x80000000L) != 0;

		prodB = prodB + (bc << 32);
		prodA = prodA + (bc >>> 32);

		bitsum = (prodB & 0x8000000000000000L) != 0;

		/* Carry */
		if ((bita && bitb) || ((bita || bitb) && !bitsum)) {
			prodA++;
		}

		/* Multiply input by shpiB */
		c = shpiB >>> 32;
		d = shpiB & 0xffffffffL;
		ac = a * c;
		bc = b * c;
		ad = a * d;

		/* Collect terms */
		ac = ac + ((bc + ad) >>> 32);

		bita = (prodB & 0x8000000000000000L) != 0;
		bitb = (ac & 0x8000000000000000L) != 0;
		prodB += ac;
		bitsum = (prodB & 0x8000000000000000L) != 0;
		/* Carry */
		if ((bita && bitb) || ((bita || bitb) && !bitsum)) {
			prodA++;
		}

		/* Multiply by shpi0 */
		c = shpi0 >>> 32;
		d = shpi0 & 0xffffffffL;

		bd = b * d;
		bc = b * c;
		ad = a * d;

		prodA += bd + ((bc + ad) << 32);

		/*
		 * prodA, prodB now contain the remainder as a fraction of PI. We want
		 * this as a fraction of PI/2, so use the following steps: 1.) multiply
		 * by 4. 2.) do a fixed point muliply by PI/4. 3.) Convert to floating
		 * point. 4.) Multiply by 2
		 */

		/* This identifies the quadrant */
		int intPart = (int) (prodA >>> 62);

		/* Multiply by 4 */
		prodA <<= 2;
		prodA |= prodB >>> 62;
		prodB <<= 2;

		/* Multiply by PI/4 */
		a = prodA >>> 32;
		b = prodA & 0xffffffffL;

		c = PI_O_4_BITS[0] >>> 32;
		d = PI_O_4_BITS[0] & 0xffffffffL;

		ac = a * c;
		bd = b * d;
		bc = b * c;
		ad = a * d;

		long prod2B = bd + (ad << 32);
		long prod2A = ac + (ad >>> 32);

		bita = (bd & 0x8000000000000000L) != 0;
		bitb = (ad & 0x80000000L) != 0;
		bitsum = (prod2B & 0x8000000000000000L) != 0;

		/* Carry */
		if ((bita && bitb) || ((bita || bitb) && !bitsum)) {
			prod2A++;
		}

		bita = (prod2B & 0x8000000000000000L) != 0;
		bitb = (bc & 0x80000000L) != 0;

		prod2B = prod2B + (bc << 32);
		prod2A = prod2A + (bc >>> 32);

		bitsum = (prod2B & 0x8000000000000000L) != 0;

		/* Carry */
		if ((bita && bitb) || ((bita || bitb) && !bitsum)) {
			prod2A++;
		}

		/* Multiply input by pio4bits[1] */
		c = PI_O_4_BITS[1] >>> 32;
		d = PI_O_4_BITS[1] & 0xffffffffL;
		ac = a * c;
		bc = b * c;
		ad = a * d;

		/* Collect terms */
		ac = ac + ((bc + ad) >>> 32);

		bita = (prod2B & 0x8000000000000000L) != 0;
		bitb = (ac & 0x8000000000000000L) != 0;
		prod2B += ac;
		bitsum = (prod2B & 0x8000000000000000L) != 0;
		/* Carry */
		if ((bita && bitb) || ((bita || bitb) && !bitsum)) {
			prod2A++;
		}

		/* Multiply inputB by pio4bits[0] */
		a = prodB >>> 32;
		b = prodB & 0xffffffffL;
		c = PI_O_4_BITS[0] >>> 32;
		d = PI_O_4_BITS[0] & 0xffffffffL;
		ac = a * c;
		bc = b * c;
		ad = a * d;

		/* Collect terms */
		ac = ac + ((bc + ad) >>> 32);

		bita = (prod2B & 0x8000000000000000L) != 0;
		bitb = (ac & 0x8000000000000000L) != 0;
		prod2B += ac;
		bitsum = (prod2B & 0x8000000000000000L) != 0;
		/* Carry */
		if ((bita && bitb) || ((bita || bitb) && !bitsum)) {
			prod2A++;
		}

		/* Convert to double */
		double tmpA = (prod2A >>> 12) / TWO_POWER_52; // High order 52 bits
		double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52; // Low
																									// bits

		double sumA = tmpA + tmpB;
		double sumB = -(sumA - tmpA - tmpB);

		/* Multiply by PI/2 and return */
		result[0] = intPart;
		result[1] = sumA * 2.0;
		result[2] = sumB * 2.0;
	}

	/**
	 * Sine function.
	 * 
	 * @param x
	 *            a number
	 * @return sin(x)
	 */
	public static double sin(double x) {
		boolean negative = false;
		int quadrant = 0;
		double xa;
		double xb = 0.0;

		/* Take absolute value of the input */
		xa = x;
		if (x < 0) {
			negative = true;
			xa = -xa;
		}

		/* Check for zero and negative zero */
		if (xa == 0.0) {
			long bits = Double.doubleToLongBits(x);
			if (bits < 0) {
				return -0.0;
			}
			return 0.0;
		}

		if (xa != xa || xa == Double.POSITIVE_INFINITY) {
			return Double.NaN;
		}

		/* Perform any argument reduction */
		if (xa > 3294198.0) {
			// PI * (2**20)
			// Argument too big for CodyWaite reduction. Must use
			// PayneHanek.
			double reduceResults[] = new double[3];
			reducePayneHanek(xa, reduceResults);
			quadrant = ((int) reduceResults[0]) & 3;
			xa = reduceResults[1];
			xb = reduceResults[2];
		} else if (xa > 1.5707963267948966) {
			/* Inline the Cody/Waite reduction for performance */

			// Estimate k
			// k = (int)(xa / 1.5707963267948966);
			int k = (int) (xa * 0.6366197723675814);

			// Compute remainder
			double remA;
			double remB;
			while (true) {
				double a = -k * 1.570796251296997;
				remA = xa + a;
				remB = -(remA - xa - a);

				a = -k * 7.549789948768648E-8;
				double b = remA;
				remA = a + b;
				remB += -(remA - b - a);

				a = -k * 6.123233995736766E-17;
				b = remA;
				remA = a + b;
				remB += -(remA - b - a);

				if (remA > 0.0)
					break;

				// Remainder is negative, so decrement k and try again.
				// This should only happen if the input is very close
				// to an even multiple of pi/2
				k--;
			}
			quadrant = k & 3;
			xa = remA;
			xb = remB;
		}

		if (negative) {
			quadrant ^= 2; // Flip bit 1
		}

		switch (quadrant) {
		case 0:
			return sinQ(xa, xb);
		case 1:
			return cosQ(xa, xb);
		case 2:
			return -sinQ(xa, xb);
		case 3:
			return -cosQ(xa, xb);
		default:
			return Double.NaN;
		}
	}

	/**
	 * Cosine function
	 * 
	 * @param x
	 *            a number
	 * @return cos(x)
	 */
	public static double cos(double x) {
		int quadrant = 0;

		/* Take absolute value of the input */
		double xa = x;
		if (x < 0) {
			xa = -xa;
		}

		if (xa != xa || xa == Double.POSITIVE_INFINITY) {
			return Double.NaN;
		}

		/* Perform any argument reduction */
		double xb = 0;
		if (xa > 3294198.0) {
			// PI * (2**20)
			// Argument too big for CodyWaite reduction. Must use
			// PayneHanek.
			double reduceResults[] = new double[3];
			reducePayneHanek(xa, reduceResults);
			quadrant = ((int) reduceResults[0]) & 3;
			xa = reduceResults[1];
			xb = reduceResults[2];
		} else if (xa > 1.5707963267948966) {
			/* Inline the Cody/Waite reduction for performance */

			// Estimate k
			// k = (int)(xa / 1.5707963267948966);
			int k = (int) (xa * 0.6366197723675814);

			// Compute remainder
			double remA;
			double remB;
			while (true) {
				double a = -k * 1.570796251296997;
				remA = xa + a;
				remB = -(remA - xa - a);

				a = -k * 7.549789948768648E-8;
				double b = remA;
				remA = a + b;
				remB += -(remA - b - a);

				a = -k * 6.123233995736766E-17;
				b = remA;
				remA = a + b;
				remB += -(remA - b - a);

				if (remA > 0.0)
					break;

				// Remainder is negative, so decrement k and try again.
				// This should only happen if the input is very close
				// to an even multiple of pi/2
				k--;
			}
			quadrant = k & 3;
			xa = remA;
			xb = remB;
		}

		// if (negative)
		// quadrant = (quadrant + 2) % 4;

		switch (quadrant) {
		case 0:
			return cosQ(xa, xb);
		case 1:
			return -sinQ(xa, xb);
		case 2:
			return -cosQ(xa, xb);
		case 3:
			return sinQ(xa, xb);
		default:
			return Double.NaN;
		}
	}

	/**
	 * Tangent function
	 * 
	 * @param x
	 *            a number
	 * @return tan(x)
	 */
	public static double tan(double x) {
		boolean negative = false;
		int quadrant = 0;

		/* Take absolute value of the input */
		double xa = x;
		if (x < 0) {
			negative = true;
			xa = -xa;
		}

		/* Check for zero and negative zero */
		if (xa == 0.0) {
			long bits = Double.doubleToLongBits(x);
			if (bits < 0) {
				return -0.0;
			}
			return 0.0;
		}

		if (xa != xa || xa == Double.POSITIVE_INFINITY) {
			return Double.NaN;
		}

		/* Perform any argument reduction */
		double xb = 0;
		if (xa > 3294198.0) {
			// PI * (2**20)
			// Argument too big for CodyWaite reduction. Must use
			// PayneHanek.
			double reduceResults[] = new double[3];
			reducePayneHanek(xa, reduceResults);
			quadrant = ((int) reduceResults[0]) & 3;
			xa = reduceResults[1];
			xb = reduceResults[2];
		} else if (xa > 1.5707963267948966) {
			/* Inline the Cody/Waite reduction for performance */

			// Estimate k
			// k = (int)(xa / 1.5707963267948966);
			int k = (int) (xa * 0.6366197723675814);

			// Compute remainder
			double remA;
			double remB;
			while (true) {
				double a = -k * 1.570796251296997;
				remA = xa + a;
				remB = -(remA - xa - a);

				a = -k * 7.549789948768648E-8;
				double b = remA;
				remA = a + b;
				remB += -(remA - b - a);

				a = -k * 6.123233995736766E-17;
				b = remA;
				remA = a + b;
				remB += -(remA - b - a);

				if (remA > 0.0)
					break;

				// Remainder is negative, so decrement k and try again.
				// This should only happen if the input is very close
				// to an even multiple of pi/2
				k--;
			}
			quadrant = k & 3;
			xa = remA;
			xb = remB;
		}

		if (xa > 1.5) {
			// Accurracy suffers between 1.5 and PI/2
			final double pi2a = 1.5707963267948966;
			final double pi2b = 6.123233995736766E-17;

			final double a = pi2a - xa;
			double b = -(a - pi2a + xa);
			b += pi2b - xb;

			xa = a + b;
			xb = -(xa - a - b);
			quadrant ^= 1;
			negative ^= true;
		}

		double result;
		if ((quadrant & 1) == 0) {
			result = tanQ(xa, xb, false);
		} else {
			result = -tanQ(xa, xb, true);
		}

		if (negative) {
			result = -result;
		}

		return result;
	}

	/**
	 * Arctangent function
	 * 
	 * @param x
	 *            a number
	 * @return atan(x)
	 */
	public static double atan(double x) {
		return atan(x, 0.0, false);
	}

	/**
	 * Internal helper function to compute arctangent.
	 * 
	 * @param xa
	 *            number from which arctangent is requested
	 * @param xb
	 *            extra bits for x (may be 0.0)
	 * @param leftPlane
	 *            if true, result angle must be put in the left half plane
	 * @return atan(xa + xb) (or angle shifted by {@code PI} if leftPlane is
	 *         true)
	 */
	private static double atan(double xa, double xb, boolean leftPlane) {
		boolean negate = false;
		int idx;

		if (xa == 0.0) { // Matches +/- 0.0; return correct sign
			return leftPlane ? copySign(Math.PI, xa) : xa;
		}

		if (xa < 0) {
			// negative
			xa = -xa;
			xb = -xb;
			negate = true;
		}

		if (xa > 1.633123935319537E16) { // Very large input
			return (negate ^ leftPlane) ? (-Math.PI / 2.0) : (Math.PI / 2.0);
		}

		/*
		 * Estimate the closest tabulated arctan value, compute eps =
		 * xa-tangentTable
		 */
		if (xa < 1.0) {
			idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5);
		} else {
			double temp = 1.0 / xa;
			idx = (int) (-((-1.7168146928204136 * temp * temp + 8.0) * temp) + 13.07);
		}
		double epsA = xa - TANGENT_TABLE_A[idx];
		double epsB = -(epsA - xa + TANGENT_TABLE_A[idx]);
		epsB += xb - TANGENT_TABLE_B[idx];

		double temp = epsA + epsB;
		epsB = -(temp - epsA - epsB);
		epsA = temp;

		/* Compute eps = eps / (1.0 + xa*tangent) */
		temp = xa * HEX_40000000;
		double ya = xa + temp - temp;
		double yb = xb + xa - ya;
		xa = ya;
		xb += yb;

		// if (idx > 8 || idx == 0)
		if (idx == 0) {
			/*
			 * If the slope of the arctan is gentle enough (< 0.45), this
			 * approximation will suffice
			 */
			// double denom = 1.0 / (1.0 + xa*tangentTableA[idx] +
			// xb*tangentTableA[idx] + xa*tangentTableB[idx] +
			// xb*tangentTableB[idx]);
			double denom = 1.0 / (1.0 + (xa + xb) * (TANGENT_TABLE_A[idx] + TANGENT_TABLE_B[idx]));
			// double denom = 1.0 / (1.0 + xa*tangentTableA[idx]);
			ya = epsA * denom;
			yb = epsB * denom;
		} else {
			double temp2 = xa * TANGENT_TABLE_A[idx];
			double za = 1.0 + temp2;
			double zb = -(za - 1.0 - temp2);
			temp2 = xb * TANGENT_TABLE_A[idx] + xa * TANGENT_TABLE_B[idx];
			temp = za + temp2;
			zb += -(temp - za - temp2);
			za = temp;

			zb += xb * TANGENT_TABLE_B[idx];
			ya = epsA / za;

			temp = ya * HEX_40000000;
			final double yaa = (ya + temp) - temp;
			final double yab = ya - yaa;

			temp = za * HEX_40000000;
			final double zaa = (za + temp) - temp;
			final double zab = za - zaa;

			/* Correct for rounding in division */
			yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za;

			yb += -epsA * zb / za / za;
			yb += epsB / za;
		}

		epsA = ya;
		epsB = yb;

		/* Evaluate polynomial */
		double epsA2 = epsA * epsA;

		/*
		 * yb = -0.09001346640161823; yb = yb * epsA2 + 0.11110718400605211; yb
		 * = yb * epsA2 + -0.1428571349122913; yb = yb * epsA2 +
		 * 0.19999999999273194; yb = yb * epsA2 + -0.33333333333333093; yb = yb
		 * * epsA2 * epsA;
		 */

		yb = 0.07490822288864472;
		yb = yb * epsA2 + -0.09088450866185192;
		yb = yb * epsA2 + 0.11111095942313305;
		yb = yb * epsA2 + -0.1428571423679182;
		yb = yb * epsA2 + 0.19999999999923582;
		yb = yb * epsA2 + -0.33333333333333287;
		yb = yb * epsA2 * epsA;

		ya = epsA;

		temp = ya + yb;
		yb = -(temp - ya - yb);
		ya = temp;

		/* Add in effect of epsB. atan'(x) = 1/(1+x^2) */
		yb += epsB / (1.0 + epsA * epsA);

		double result;
		double resultb;

		// result = yb + eighths[idx] + ya;
		double za = EIGHTHS[idx] + ya;
		double zb = -(za - EIGHTHS[idx] - ya);
		temp = za + yb;
		zb += -(temp - za - yb);
		za = temp;

		result = za + zb;
		resultb = -(result - za - zb);

		if (leftPlane) {
			// Result is in the left plane
			final double pia = 1.5707963267948966 * 2.0;
			final double pib = 6.123233995736766E-17 * 2.0;

			za = pia - result;
			zb = -(za - pia + result);
			zb += pib - resultb;

			result = za + zb;
			resultb = -(result - za - zb);
		}

		if (negate ^ leftPlane) {
			result = -result;
		}

		return result;
	}

	/**
	 * Two arguments arctangent function
	 * 
	 * @param y
	 *            ordinate
	 * @param x
	 *            abscissa
	 * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
	 */
	public static double atan2(double y, double x) {
		if (x != x || y != y) {
			return Double.NaN;
		}

		if (y == 0.0) {
			double result = x * y;
			double invx = 1.0 / x;
			double invy = 1.0 / y;

			if (invx == 0.0) { // X is infinite
				if (x > 0) {
					return y; // return +/- 0.0
				} else {
					return copySign(Math.PI, y);
				}
			}

			if (x < 0.0 || invx < 0.0) {
				if (y < 0.0 || invy < 0.0) {
					return -Math.PI;
				} else {
					return Math.PI;
				}
			} else {
				return result;
			}
		}

		// y cannot now be zero

		if (y == Double.POSITIVE_INFINITY) {
			if (x == Double.POSITIVE_INFINITY) {
				return Math.PI / 4.0;
			}

			if (x == Double.NEGATIVE_INFINITY) {
				return Math.PI * 3.0 / 4.0;
			}

			return Math.PI / 2.0;
		}

		if (y == Double.NEGATIVE_INFINITY) {
			if (x == Double.POSITIVE_INFINITY) {
				return -Math.PI / 4.0;
			}

			if (x == Double.NEGATIVE_INFINITY) {
				return -Math.PI * 3.0 / 4.0;
			}

			return -Math.PI / 2.0;
		}

		if (x == Double.POSITIVE_INFINITY) {
			if (y > 0.0 || 1 / y > 0.0) {
				return 0.0;
			}

			if (y < 0.0 || 1 / y < 0.0) {
				return -0.0;
			}
		}

		if (x == Double.NEGATIVE_INFINITY) {
			if (y > 0.0 || 1 / y > 0.0) {
				return Math.PI;
			}

			if (y < 0.0 || 1 / y < 0.0) {
				return -Math.PI;
			}
		}

		// Neither y nor x can be infinite or NAN here

		if (x == 0) {
			if (y > 0.0 || 1 / y > 0.0) {
				return Math.PI / 2.0;
			}

			if (y < 0.0 || 1 / y < 0.0) {
				return -Math.PI / 2.0;
			}
		}

		// Compute ratio r = y/x
		final double r = y / x;
		if (Double.isInfinite(r)) { // bypass calculations that can create NaN
			return atan(r, 0, x < 0);
		}

		double ra = doubleHighPart(r);
		double rb = r - ra;

		// Split x
		final double xa = doubleHighPart(x);
		final double xb = x - xa;

		rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;

		double temp = ra + rb;
		rb = -(temp - ra - rb);
		ra = temp;

		if (ra == 0) { // Fix up the sign so atan works correctly
			ra = copySign(0.0, y);
		}

		// Call atan
		double result = atan(ra, rb, x < 0);

		return result;
	}

	/**
	 * Compute the arc sine of a number.
	 * 
	 * @param x
	 *            number on which evaluation is done
	 * @return arc sine of x
	 */
	public static double asin(double x) {
		if (x != x) {
			return Double.NaN;
		}

		if (x > 1.0 || x < -1.0) {
			return Double.NaN;
		}

		if (x == 1.0) {
			return Math.PI / 2.0;
		}

		if (x == -1.0) {
			return -Math.PI / 2.0;
		}

		if (x == 0.0) { // Matches +/- 0.0; return correct sign
			return x;
		}

		/* Compute asin(x) = atan(x/sqrt(1-x*x)) */

		/* Split x */
		double temp = x * HEX_40000000;
		final double xa = x + temp - temp;
		final double xb = x - xa;

		/* Square it */
		double ya = xa * xa;
		double yb = xa * xb * 2.0 + xb * xb;

		/* Subtract from 1 */
		ya = -ya;
		yb = -yb;

		double za = 1.0 + ya;
		double zb = -(za - 1.0 - ya);

		temp = za + yb;
		zb += -(temp - za - yb);
		za = temp;

		/* Square root */
		double y;
		y = sqrt(za);
		temp = y * HEX_40000000;
		ya = y + temp - temp;
		yb = y - ya;

		/* Extend precision of sqrt */
		yb += (za - ya * ya - 2 * ya * yb - yb * yb) / (2.0 * y);

		/* Contribution of zb to sqrt */
		double dx = zb / (2.0 * y);

		// Compute ratio r = x/y
		double r = x / y;
		temp = r * HEX_40000000;
		double ra = r + temp - temp;
		double rb = r - ra;

		rb += (x - ra * ya - ra * yb - rb * ya - rb * yb) / y; // Correct for
																// rounding in
																// division
		rb += -x * dx / y / y; // Add in effect additional bits of sqrt.

		temp = ra + rb;
		rb = -(temp - ra - rb);
		ra = temp;

		return atan(ra, rb, false);
	}

	/**
	 * Compute the arc cosine of a number.
	 * 
	 * @param x  number on which evaluation is done
	 * @return arc cosine of x
	 */
	public static double acos(double x) {
		if (x != x) {
			return Double.NaN;
		} else if (x >= 1) {
			return 0;
		} else if (x <= -1) {
			return Math.PI;
		} else if (x == 0) {
			return Math.PI / 2.0;
		}


		/* Compute acos(x) = atan(sqrt(1-x*x)/x) */

		/* Split x */
		double temp = x * HEX_40000000;
		final double xa = x + temp - temp;
		final double xb = x - xa;

		/* Square it */
		double ya = xa * xa;
		double yb = xa * xb * 2.0 + xb * xb;

		/* Subtract from 1 */
		ya = -ya;
		yb = -yb;

		double za = 1.0 + ya;
		double zb = -(za - 1.0 - ya);

		temp = za + yb;
		zb += -(temp - za - yb);
		za = temp;

		/* Square root */
		double y = sqrt(za);
		temp = y * HEX_40000000;
		ya = y + temp - temp;
		yb = y - ya;

		/* Extend precision of sqrt */
		yb += (za - ya * ya - 2 * ya * yb - yb * yb) / (2.0 * y);

		/* Contribution of zb to sqrt */
		yb += zb / (2.0 * y);
		y = ya + yb;
		yb = -(y - ya - yb);

		// Compute ratio r = y/x
		double r = y / x;

		// Did r overflow?
		if (Double.isInfinite(r)) { // x is effectively zero
			return Math.PI / 2; // so return the appropriate value
		}

		double ra = doubleHighPart(r);
		double rb = r - ra;

		rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x; // Correct for
																// rounding in
																// division
		rb += yb / x; // Add in effect additional bits of sqrt.

		temp = ra + rb;
		rb = -(temp - ra - rb);
		ra = temp;

		return atan(ra, rb, x < 0);
	}

	/**
	 * Compute the cubic root of a number.
	 * 
	 * @param x
	 *            number on which evaluation is done
	 * @return cubic root of x
	 */
	public static double cbrt(double x) {
		/* Convert input double to bits */
		long inbits = Double.doubleToLongBits(x);
		int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
		boolean subnormal = false;

		if (exponent == -1023) {
			if (x == 0) {
				return x;
			}

			/* Subnormal, so normalize */
			subnormal = true;
			x *= 1.8014398509481984E16; // 2^54
			inbits = Double.doubleToLongBits(x);
			exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
		}

		if (exponent == 1024) {
			// Nan or infinity. Don't care which.
			return x;
		}

		/* Divide the exponent by 3 */
		int exp3 = exponent / 3;

		/* p2 will be the nearest power of 2 to x with its exponent divided by 3 */
		double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) | (long) (((exp3 + 1023) & 0x7ff)) << 52);

		/* This will be a number between 1 and 2 */
		final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);

		/* Estimate the cube root of mant by polynomial */
		double est = -0.010714690733195933;
		est = est * mant + 0.0875862700108075;
		est = est * mant + -0.3058015757857271;
		est = est * mant + 0.7249995199969751;
		est = est * mant + 0.5039018405998233;

		est *= CBRTTWO[exponent % 3 + 2];

		// est should now be good to about 15 bits of precision. Do 2 rounds of
		// Newton's method to get closer, this should get us full double
		// precision
		// Scale down x for the purpose of doing newtons method. This avoids
		// over/under flows.
		final double xs = x / (p2 * p2 * p2);
		est += (xs - est * est * est) / (3 * est * est);
		est += (xs - est * est * est) / (3 * est * est);

		// Do one round of Newton's method in extended precision to get the last
		// bit right.
		double temp = est * HEX_40000000;
		double ya = est + temp - temp;
		double yb = est - ya;

		double za = ya * ya;
		double zb = ya * yb * 2.0 + yb * yb;
		temp = za * HEX_40000000;
		double temp2 = za + temp - temp;
		zb += za - temp2;
		za = temp2;

		zb = za * yb + ya * zb + zb * yb;
		za = za * ya;

		double na = xs - za;
		double nb = -(na - xs + za);
		nb -= zb;

		est += (na + nb) / (3 * est * est);

		/* Scale by a power of two, so this is exact. */
		est *= p2;

		if (subnormal) {
			est *= 3.814697265625E-6; // 2^-18
		}

		return est;
	}

	/**
	 * Convert degrees to radians, with error of less than 0.5 ULP
	 * 
	 * @param x
	 *            angle in degrees
	 * @return x converted into radians
	 */
	public static double toRadians(double x) {
		if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return
												// correct sign
			return x;
		}

		// These are PI/180 split into high and low order bits
		final double facta = 0.01745329052209854;
		final double factb = 1.997844754509471E-9;

		double xa = doubleHighPart(x);
		double xb = x - xa;

		double result = xb * factb + xb * facta + xa * factb + xa * facta;
		if (result == 0) {
			result = result * x; // ensure correct sign if calculation
									// underflows
		}
		return result;
	}

	/**
	 * Convert radians to degrees, with error of less than 0.5 ULP
	 * 
	 * @param x
	 *            angle in radians
	 * @return x converted into degrees
	 */
	public static double toDegrees(double x) {
		if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return
												// correct sign
			return x;
		}

		// These are 180/PI split into high and low order bits
		final double facta = 57.2957763671875;
		final double factb = 3.145894820876798E-6;

		double xa = doubleHighPart(x);
		double xb = x - xa;

		return xb * factb + xb * facta + xa * factb + xa * facta;
	}

	/**
	 * Absolute value.
	 * 
	 * @param x
	 *            number from which absolute value is requested
	 * @return abs(x)
	 */
	public static int abs(final int x) {
		return (x < 0) ? -x : x;
	}

	/**
	 * Absolute value.
	 * 
	 * @param x
	 *            number from which absolute value is requested
	 * @return abs(x)
	 */
	public static long abs(final long x) {
		return (x < 0l) ? -x : x;
	}

	/**
	 * Absolute value.
	 * 
	 * @param x
	 *            number from which absolute value is requested
	 * @return abs(x)
	 */
	public static float abs(final float x) {
		return (x < 0.0f) ? -x : (x == 0.0f) ? 0.0f : x; // -0.0 => +0.0
	}

	/**
	 * Absolute value.
	 * 
	 * @param x
	 *            number from which absolute value is requested
	 * @return abs(x)
	 */
	public static double abs(double x) {
		return (x < 0.0) ? -x : (x == 0.0) ? 0.0 : x; // -0.0 => +0.0
	}

	/**
	 * Compute least significant bit (Unit in Last Position) for a number.
	 * 
	 * @param x
	 *            number from which ulp is requested
	 * @return ulp(x)
	 */
	public static double ulp(double x) {
		if (Double.isInfinite(x)) {
			return Double.POSITIVE_INFINITY;
		}
		return abs(x - Double.longBitsToDouble(Double.doubleToLongBits(x) ^ 1));
	}

	/**
	 * Compute least significant bit (Unit in Last Position) for a number.
	 * 
	 * @param x
	 *            number from which ulp is requested
	 * @return ulp(x)
	 */
	public static float ulp(float x) {
		if (Float.isInfinite(x)) {
			return Float.POSITIVE_INFINITY;
		}
		return abs(x - Float.intBitsToFloat(Float.floatToIntBits(x) ^ 1));
	}

	/**
	 * Multiply a double number by a power of 2.
	 * 
	 * @param d
	 *            number to multiply
	 * @param n
	 *            power of 2
	 * @return d &times; 2<sup>n</sup>
	 */
	public static double scalb(final double d, final int n) {

		// first simple and fast handling when 2^n can be represented using
		// normal numbers
		if ((n > -1023) && (n < 1024)) {
			return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
		}

		// handle special cases
		if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) {
			return d;
		}
		if (n < -2098) {
			return (d > 0) ? 0.0 : -0.0;
		}
		if (n > 2097) {
			return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
		}

		// decompose d
		final long bits = Double.doubleToLongBits(d);
		final long sign = bits & 0x8000000000000000L;
		int exponent = ((int) (bits >>> 52)) & 0x7ff;
		long mantissa = bits & 0x000fffffffffffffL;

		// compute scaled exponent
		int scaledExponent = exponent + n;

		if (n < 0) {
			// we are really in the case n <= -1023
			if (scaledExponent > 0) {
				// both the input and the result are normal numbers, we only
				// adjust the exponent
				return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
			} else if (scaledExponent > -53) {
				// the input is a normal number and the result is a subnormal
				// number

				// recover the hidden mantissa bit
				mantissa = mantissa | (1L << 52);

				// scales down complete mantissa, hence losing least significant
				// bits
				final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
				mantissa = mantissa >>> (1 - scaledExponent);
				if (mostSignificantLostBit != 0) {
					// we need to add 1 bit to round up the result
					mantissa++;
				}
				return Double.longBitsToDouble(sign | mantissa);

			} else {
				// no need to compute the mantissa, the number scales down to 0
				return (sign == 0L) ? 0.0 : -0.0;
			}
		} else {
			// we are really in the case n >= 1024
			if (exponent == 0) {

				// the input number is subnormal, normalize it
				while ((mantissa >>> 52) != 1) {
					mantissa = mantissa << 1;
					--scaledExponent;
				}
				++scaledExponent;
				mantissa = mantissa & 0x000fffffffffffffL;

				if (scaledExponent < 2047) {
					return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
				} else {
					return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
				}

			} else if (scaledExponent < 2047) {
				return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
			} else {
				return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
			}
		}

	}

	/**
	 * Multiply a float number by a power of 2.
	 * 
	 * @param f
	 *            number to multiply
	 * @param n
	 *            power of 2
	 * @return f &times; 2<sup>n</sup>
	 */
	public static float scalb(final float f, final int n) {

		// first simple and fast handling when 2^n can be represented using
		// normal numbers
		if ((n > -127) && (n < 128)) {
			return f * Float.intBitsToFloat((n + 127) << 23);
		}

		// handle special cases
		if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) {
			return f;
		}
		if (n < -277) {
			return (f > 0) ? 0.0f : -0.0f;
		}
		if (n > 276) {
			return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
		}

		// decompose f
		final int bits = Float.floatToIntBits(f);
		final int sign = bits & 0x80000000;
		int exponent = (bits >>> 23) & 0xff;
		int mantissa = bits & 0x007fffff;

		// compute scaled exponent
		int scaledExponent = exponent + n;

		if (n < 0) {
			// we are really in the case n <= -127
			if (scaledExponent > 0) {
				// both the input and the result are normal numbers, we only
				// adjust the exponent
				return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
			} else if (scaledExponent > -24) {
				// the input is a normal number and the result is a subnormal
				// number

				// recover the hidden mantissa bit
				mantissa = mantissa | (1 << 23);

				// scales down complete mantissa, hence losing least significant
				// bits
				final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent));
				mantissa = mantissa >>> (1 - scaledExponent);
				if (mostSignificantLostBit != 0) {
					// we need to add 1 bit to round up the result
					mantissa++;
				}
				return Float.intBitsToFloat(sign | mantissa);

			} else {
				// no need to compute the mantissa, the number scales down to 0
				return (sign == 0) ? 0.0f : -0.0f;
			}
		} else {
			// we are really in the case n >= 128
			if (exponent == 0) {

				// the input number is subnormal, normalize it
				while ((mantissa >>> 23) != 1) {
					mantissa = mantissa << 1;
					--scaledExponent;
				}
				++scaledExponent;
				mantissa = mantissa & 0x007fffff;

				if (scaledExponent < 255) {
					return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
				} else {
					return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
				}

			} else if (scaledExponent < 255) {
				return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
			} else {
				return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
			}
		}

	}

	/**
	 * Get the next machine representable number after a number, moving in the
	 * direction of another number.
	 * <p>
	 * The ordering is as follows (increasing):
	 * <ul>
	 * <li>-INFINITY</li>
	 * <li>-MAX_VALUE</li>
	 * <li>-MIN_VALUE</li>
	 * <li>-0.0</li>
	 * <li>+0.0</li>
	 * <li>+MIN_VALUE</li>
	 * <li>+MAX_VALUE</li>
	 * <li>+INFINITY</li>
	 * <li></li>
	 * <p>
	 * If arguments compare equal, then the second argument is returned.
	 * <p>
	 * If {@code direction} is greater than {@code d}, the smallest machine
	 * representable number strictly greater than {@code d} is returned; if
	 * less, then the largest representable number strictly less than {@code d}
	 * is returned.
	 * </p>
	 * <p>
	 * If {@code d} is infinite and direction does not bring it back to finite
	 * numbers, it is returned unchanged.
	 * </p>
	 *
	 * @param d
	 *            base number
	 * @param direction
	 *            (the only important thing is whether {@code direction} is
	 *            greater or smaller than {@code d})
	 * @return the next machine representable number in the specified direction
	 */
	public static double nextAfter(double d, double direction) {

		// handling of some important special cases
		if (Double.isNaN(d) || Double.isNaN(direction)) {
			return Double.NaN;
		} else if (d == direction) {
			return direction;
		} else if (Double.isInfinite(d)) {
			return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE;
		} else if (d == 0) {
			return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
		}
		// special cases MAX_VALUE to infinity and MIN_VALUE to 0
		// are handled just as normal numbers

		final long bits = Double.doubleToLongBits(d);
		final long sign = bits & 0x8000000000000000L;
		if ((direction < d) ^ (sign == 0L)) {
			return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1));
		} else {
			return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1));
		}

	}

	/**
	 * Get the next machine representable number after a number, moving in the
	 * direction of another number.
	 * <p>
	 * The ordering is as follows (increasing):
	 * <ul>
	 * <li>-INFINITY</li>
	 * <li>-MAX_VALUE</li>
	 * <li>-MIN_VALUE</li>
	 * <li>-0.0</li>
	 * <li>+0.0</li>
	 * <li>+MIN_VALUE</li>
	 * <li>+MAX_VALUE</li>
	 * <li>+INFINITY</li>
	 * <li></li>
	 * <p>
	 * If arguments compare equal, then the second argument is returned.
	 * <p>
	 * If {@code direction} is greater than {@code f}, the smallest machine
	 * representable number strictly greater than {@code f} is returned; if
	 * less, then the largest representable number strictly less than {@code f}
	 * is returned.
	 * </p>
	 * <p>
	 * If {@code f} is infinite and direction does not bring it back to finite
	 * numbers, it is returned unchanged.
	 * </p>
	 *
	 * @param f
	 *            base number
	 * @param direction
	 *            (the only important thing is whether {@code direction} is
	 *            greater or smaller than {@code f})
	 * @return the next machine representable number in the specified direction
	 */
	public static float nextAfter(final float f, final double direction) {

		// handling of some important special cases
		if (Double.isNaN(f) || Double.isNaN(direction)) {
			return Float.NaN;
		} else if (f == direction) {
			return (float) direction;
		} else if (Float.isInfinite(f)) {
			return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE;
		} else if (f == 0f) {
			return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE;
		}
		// special cases MAX_VALUE to infinity and MIN_VALUE to 0
		// are handled just as normal numbers

		final int bits = Float.floatToIntBits(f);
		final int sign = bits & 0x80000000;
		if ((direction < f) ^ (sign == 0)) {
			return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1));
		} else {
			return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1));
		}

	}

	/**
	 * Get the largest whole number smaller than x.
	 * 
	 * @param x
	 *            number from which floor is requested
	 * @return a double number f such that f is an integer f <= x < f + 1.0
	 */
	public static double floor(double x) {
		long y;

		if (x != x) { // NaN
			return x;
		}

		if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) {
			return x;
		}

		y = (long) x;
		if (x < 0 && y != x) {
			y--;
		}

		if (y == 0) {
			return x * y;
		}

		return y;
	}

	/**
	 * Get the smallest whole number larger than x.
	 * 
	 * @param x
	 *            number from which ceil is requested
	 * @return a double number c such that c is an integer c - 1.0 < x <= c
	 */
	public static double ceil(double x) {
		double y;

		if (x != x) { // NaN
			return x;
		}

		y = floor(x);
		if (y == x) {
			return y;
		}

		y += 1.0;

		if (y == 0) {
			return x * y;
		}

		return y;
	}

	/**
	 * Get the whole number that is the nearest to x, or the even one if x is
	 * exactly half way between two integers.
	 * 
	 * @param x
	 *            number from which nearest whole number is requested
	 * @return a double number r such that r is an integer r - 0.5 <= x <= r +
	 *         0.5
	 */
	public static double rint(double x) {
		double y = floor(x);
		double d = x - y;

		if (d > 0.5) {
			if (y == -1.0) {
				return -0.0; // Preserve sign of operand
			}
			return y + 1.0;
		}
		if (d < 0.5) {
			return y;
		}

		/* half way, round to even */
		long z = (long) y;
		return (z & 1) == 0 ? y : y + 1.0;
	}

	/**
	 * Get the closest long to x.
	 * 
	 * @param x
	 *            number from which closest long is requested
	 * @return closest long to x
	 */
	public static long round(double x) {
		return (long) floor(x + 0.5);
	}

	/**
	 * Get the closest int to x.
	 * 
	 * @param x
	 *            number from which closest int is requested
	 * @return closest int to x
	 */
	public static int round(final float x) {
		return (int) floor(x + 0.5f);
	}

	/**
	 * Compute the minimum of two values
	 * 
	 * @param a
	 *            first value
	 * @param b
	 *            second value
	 * @return a if a is lesser or equal to b, b otherwise
	 */
	public static int min(final int a, final int b) {
		return (a <= b) ? a : b;
	}

	/**
	 * Compute the minimum of two values
	 * 
	 * @param a
	 *            first value
	 * @param b
	 *            second value
	 * @return a if a is lesser or equal to b, b otherwise
	 */
	public static long min(final long a, final long b) {
		return (a <= b) ? a : b;
	}

	/**
	 * Compute the minimum of two values
	 * 
	 * @param a
	 *            first value
	 * @param b
	 *            second value
	 * @return a if a is lesser or equal to b, b otherwise
	 */
	public static float min(final float a, final float b) {
		if (a > b) {
			return b;
		}
		if (a < b) {
			return a;
		}
		/* if either arg is NaN, return NaN */
		if (a != b) {
			return Float.NaN;
		}
		/* min(+0.0,-0.0) == -0.0 */
		/* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
		int bits = Float.floatToRawIntBits(a);
		if (bits == 0x80000000) {
			return a;
		}
		return b;
	}

	/**
	 * Compute the minimum of two values
	 * 
	 * @param a
	 *            first value
	 * @param b
	 *            second value
	 * @return a if a is lesser or equal to b, b otherwise
	 */
	public static double min(final double a, final double b) {
		if (a > b) {
			return b;
		}
		if (a < b) {
			return a;
		}
		/* if either arg is NaN, return NaN */
		if (a != b) {
			return Double.NaN;
		}
		/* min(+0.0,-0.0) == -0.0 */
		/* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
		long bits = Double.doubleToRawLongBits(a);
		if (bits == 0x8000000000000000L) {
			return a;
		}
		return b;
	}

	/**
	 * Compute the maximum of two values
	 * 
	 * @param a
	 *            first value
	 * @param b
	 *            second value
	 * @return b if a is lesser or equal to b, a otherwise
	 */
	public static int max(final int a, final int b) {
		return (a <= b) ? b : a;
	}

	/**
	 * Compute the maximum of two values
	 * 
	 * @param a
	 *            first value
	 * @param b
	 *            second value
	 * @return b if a is lesser or equal to b, a otherwise
	 */
	public static long max(final long a, final long b) {
		return (a <= b) ? b : a;
	}

	/**
	 * Compute the maximum of two values
	 * 
	 * @param a
	 *            first value
	 * @param b
	 *            second value
	 * @return b if a is lesser or equal to b, a otherwise
	 */
	public static float max(final float a, final float b) {
		if (a > b) {
			return a;
		}
		if (a < b) {
			return b;
		}
		/* if either arg is NaN, return NaN */
		if (a != b) {
			return Float.NaN;
		}
		/* min(+0.0,-0.0) == -0.0 */
		/* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
		int bits = Float.floatToRawIntBits(a);
		if (bits == 0x80000000) {
			return b;
		}
		return a;
	}

	/**
	 * Compute the maximum of two values
	 * 
	 * @param a
	 *            first value
	 * @param b
	 *            second value
	 * @return b if a is lesser or equal to b, a otherwise
	 */
	public static double max(final double a, final double b) {
		if (a > b) {
			return a;
		}
		if (a < b) {
			return b;
		}
		/* if either arg is NaN, return NaN */
		if (a != b) {
			return Double.NaN;
		}
		/* min(+0.0,-0.0) == -0.0 */
		/* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
		long bits = Double.doubleToRawLongBits(a);
		if (bits == 0x8000000000000000L) {
			return b;
		}
		return a;
	}

	/**
	 * Returns the hypotenuse of a triangle with sides {@code x} and {@code y} -
	 * sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)<br/>
	 * avoiding intermediate overflow or underflow.
	 *
	 * <ul>
	 * <li>If either argument is infinite, then the result is positive infinity.
	 * </li>
	 * <li>else, if either argument is NaN then the result is NaN.</li>
	 * </ul>
	 *
	 * @param x
	 *            a value
	 * @param y
	 *            a value
	 * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
	 */
	public static double hypot(final double x, final double y) {
		if (Double.isInfinite(x) || Double.isInfinite(y)) {
			return Double.POSITIVE_INFINITY;
		} else if (Double.isNaN(x) || Double.isNaN(y)) {
			return Double.NaN;
		} else {

			final int expX = getExponent(x);
			final int expY = getExponent(y);
			if (expX > expY + 27) {
				// y is neglectible with respect to x
				return abs(x);
			} else if (expY > expX + 27) {
				// x is neglectible with respect to y
				return abs(y);
			} else {

				// find an intermediate scale to avoid both overflow and
				// underflow
				final int middleExp = (expX + expY) / 2;

				// scale parameters without losing precision
				final double scaledX = scalb(x, -middleExp);
				final double scaledY = scalb(y, -middleExp);

				// compute scaled hypotenuse
				final double scaledH = sqrt(scaledX * scaledX + scaledY * scaledY);

				// remove scaling
				return scalb(scaledH, middleExp);

			}

		}
	}

	/**
	 * Computes the remainder as prescribed by the IEEE 754 standard. The
	 * remainder value is mathematically equal to {@code x - y*n} where
	 * {@code n} is the mathematical integer closest to the exact mathematical
	 * value of the quotient {@code x/y}. If two mathematical integers are
	 * equally close to {@code x/y} then {@code n} is the integer that is even.
	 * <p>
	 * <ul>
	 * <li>If either operand is NaN, the result is NaN.</li>
	 * <li>If the result is not NaN, the sign of the result equals the sign of
	 * the dividend.</li>
	 * <li>If the dividend is an infinity, or the divisor is a zero, or both,
	 * the result is NaN.</li>
	 * <li>If the dividend is finite and the divisor is an infinity, the result
	 * equals the dividend.</li>
	 * <li>If the dividend is a zero and the divisor is finite, the result
	 * equals the dividend.</li>
	 * </ul>
	 * <p>
	 * <b>Note:</b> this implementation currently delegates to
	 * {@link StrictMath#IEEEremainder}
	 * 
	 * @param dividend
	 *            the number to be divided
	 * @param divisor
	 *            the number by which to divide
	 * @return the remainder, rounded
	 */
	public static double IEEEremainder(double dividend, double divisor) {
		return StrictMath.IEEEremainder(dividend, divisor); // TODO provide our
															// own
															// implementation
	}

	/**
	 * Returns the first argument with the sign of the second argument. A NaN
	 * {@code sign} argument is treated as positive.
	 *
	 * @param magnitude
	 *            the value to return
	 * @param sign
	 *            the sign for the returned value
	 * @return the magnitude with the same sign as the {@code sign} argument
	 */
	public static double copySign(double magnitude, double sign) {
		long m = Double.doubleToLongBits(magnitude);
		long s = Double.doubleToLongBits(sign);
		if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
			return magnitude;
		}
		return -magnitude; // flip sign
	}

	/**
	 * Returns the first argument with the sign of the second argument. A NaN
	 * {@code sign} argument is treated as positive.
	 *
	 * @param magnitude
	 *            the value to return
	 * @param sign
	 *            the sign for the returned value
	 * @return the magnitude with the same sign as the {@code sign} argument
	 */
	public static float copySign(float magnitude, float sign) {
		int m = Float.floatToIntBits(magnitude);
		int s = Float.floatToIntBits(sign);
		if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
			return magnitude;
		}
		return -magnitude; // flip sign
	}

	/**
	 * Return the exponent of a double number, removing the bias.
	 * <p>
	 * For double numbers of the form 2<sup>x</sup>, the unbiased exponent is
	 * exactly x.
	 * </p>
	 * 
	 * @param d
	 *            number from which exponent is requested
	 * @return exponent for d in IEEE754 representation, without bias
	 */
	public static int getExponent(final double d) {
		return (int) ((Double.doubleToLongBits(d) >>> 52) & 0x7ff) - 1023;
	}

	/**
	 * Return the exponent of a float number, removing the bias.
	 * <p>
	 * For float numbers of the form 2<sup>x</sup>, the unbiased exponent is
	 * exactly x.
	 * </p>
	 * 
	 * @param f
	 *            number from which exponent is requested
	 * @return exponent for d in IEEE754 representation, without bias
	 */
	public static int getExponent(final float f) {
		return ((Float.floatToIntBits(f) >>> 23) & 0xff) - 127;
	}

	public static void main(String[] args) {
		for (int k = 0; k < 10; k++) {
			{
				long s = System.currentTimeMillis();
				double sum = 0;
				for (double i = -1; i < 1; i += .000001) {
					sum += Math.acos(i);
				}
				System.out.println("1Sum = " + sum);
				System.out.println("1Done= " + (System.currentTimeMillis() - s));
			}
			{
				long s = System.currentTimeMillis();
				double sum = 0;
				for (double i = -1; i < 1; i += .000001) {
					sum += acos(i);
				}
				System.out.println("2Sum = " + sum);
				System.out.println("2Done= " + (System.currentTimeMillis() - s));
			}
		}

	}
}
